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Abstract: 

The treatment of Heisenberg films with many-body Green's function theory (GFT) 
is reviewed. The basic equations of GFT are derived in sufficient detail so that the 
rest of the paper can be understood without having to consult further literature. 
The main part of the paper is concerned with applications of the formalism to fer- 
romagnetic, antiferromagnetic and coupled ferromagnetic-antiferromagnetic Heisen- 
berg films based on a generalized Tyablikov (RPA) decoupling of the exchange in- 
teraction and exchange anisotropy terms and an Anderson-Callen decoupling for a 
weak single-ion anisotropy. We not only give a consistent description of our own 
work but also refer extensively to related investigations. We discuss in particular 
the reorientation of the magnetization as a function of the temperature and film 
thickness. If the single-ion anisotropy is strong, it can be treated exactly by go- 
ing to higher-order Green's functions. We also discuss the extension of the theory 
beyond RPA. Finally the limitations of GFT are pointed out. 
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1. Introduction and outline 

Many-body GF theory is used in many fields in statistical mechanics ( e.g. see the 
early reviews [T| or j2j, or the more recent book [3 ). Extensive applications of the 
formalism to the theory of magnetism can be found in the books |lj and G . In 
the sixties and seventies of the last century, emphasis was put on the properties of 
bulk magnets. Since then, the advance in experimental techniques has stimulated 
an increasing interest in magnetic systems with reduced dimension. 

One main stream in current research is the attempt to describe 3D magnetic 
systems with strong electron-electron correlations in terms of electronic structures 
with the help of ab initio calculations: Density Functional Theory (DFT), which 
is successful by itself for systems with weak electron-electron correlations, must 
here be combined with many-body techniques, as in Dynamical Mean Field Theory 
(DMFT). For a recent review see e.g. Ref. jH]. 

The present paper is concerned with less ambitious models based on the Heisen- 
berg Hamiltonian with the inclusion of anisotropics. It provides an overview of 
many-body Green's function (GF) techniques applied to the magnetic properties of 
layered two-dimensional structures; i.e. it is concerned essentially with quasi-two- 
dimensional Heisenberg films. The techniques developed in the present paper may 
also be useful for treating cases in which a Heisenberg kind of Hamiltonian is derived 
from a microscopic theory. Emphasis is put on the development of the formalism. 
For a paper that discusses the relevant experimental situation in more detail we 
refer to Ref. [7j. 

In Section 2, we derive the direct Heisenberg exchange interaction by plau- 
sibility arguments: an orthogonal basis leads to ferromagnetic exchange and a 
non-orthogonal basis in the Heitler-London framework allows antiferromagnetic ex- 
change. In Section 3, we derive the basic equations of the formalism for the double- 
time Green's function in sufficient detail that it should not be necessary to con- 
sult any further literature to understand the rest of the paper. Section 4 deals 
with applications of the GF formalism to Heisenberg films. In the pedagogical Sec- 
tion 4.1 a ferromagnetic spin 5 = 1/2 Heisenberg monolayer in a magnetic field is 
treated. Section 4.2 deals with ferromagnetic Heisenberg films with anisotropies for 
general spin and Section 4.3 considers antiferromagnetic and coupled ferromagnetic- 
antiferromagnetic Heisenberg films. Section 5 extends the formalism beyond the 
RPA approach. Section 6 presents our conclusions and points out some open prob- 
lems and limitations of Green's function theory (GFT). 
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2. The Heisenberg exchange interaction 

The present article describes magnetic systems in terms of a Heisenberg Hamiltonian 
and anisotropy terms. This is a phenomenological approach, in which the strengths 
of the exchange interaction and anisotropies are considered as parameters which 
could be fitted to experiments. In this section, we discuss the origin of the Heisenberg 
exchange Hamiltonian. 

The exchange interaction is a manifestation of the Coulomb interaction and 
quantum-mechanical indistinguishability (the Pauli principle). It is quite compli- 
cated to derive exchange Hamiltonians from first principles, but this is often pos- 
sible by adopting adequate approximations. The form of the Heisenberg exchange 
Hamiltonian is a gross simplification that can, however, be made plausible for simple 
cases. The ferromagnetic direct exchange can be derived from a two-electron model 
by assuming orthogonal basis states. The direct antiferromagnetic exchange can be 
made plausible with the Heitler-London scheme. 

2.1. Direct exchange with orthogonal basis states (ferromagnetism) 

Consider two electrons (e.g. in a 3d 2 configuration) and a Hamiltonian consisting 
of the sum of two single-electron Hamiltonians, /io( r i) and foo( r 2) and the Coulomb 
interaction: 

T~t = ^0(1*1) + h (r 2 ) + 1 — r, (1) 

— r2 1 

where the single-electron problem is assumed to be solved: 

h (r)(f) a (r) = e a (f) a (r), 
(0 a (r)|0 6 (r)> = 5 ab . (2) 

The electrons can couple to triplet (S — 1) or singlet (S = 0) states with wave 
functions characterized by \S, S z ): 

|1,1> =?/>!, |1,0) = i=(^ 2 + ^3), |1,-1)=^4, 

10,0) = -L(^-^s), (3) 

and 



Vl = 


^Xt( s i)Xt(s2)[0a(ri)06(r 2 ) 


- 0a(r 2 )0 b (ri)] , 




^2 = 


^[X|(si)x T ( s 2)0a(ri)<Mr 2 ) 


- Xt( s i)X|(s2)<Mri)0 a (r 2 )] , 




^3 = 


^[Xt(si)X|(s2)0a(ri)(Mr 2 ) 


- X|(si)Xt( s 2)<Mri)0„(r 2 )] , 




^4 = 


^X|(si)X|(s 2 )[0a(ri)06(r 2 ) 


- 0a(r 2 )0 b (r 1 )] . 


(4) 
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Here Xui) are the spin wave functions with spins up or down. 
Defining a Coulomb integal as 



C ab = e 2 J dr 1 j 
and an exchange integral as 

Jab = e 2 J dri J 



dr. 



|0q(ri)| 2 |0 b (r 2 )| 
ki - r 2 | 



dr? 



i*(ri)&(ri)$;(r 2 )0 o (r 2 ) 



|ri - r 2 | 

we may express the matrix elements of the Coulomb interaction as 

C a b Jab j 
\lfa) = Cab , 
— —Jab ■ 
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|ri - r 2 
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|ri - r 2 
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Jab J 



We then find for the Hamiltonian matrix 



(e + Q>) 



There are three degenerate eigenvalues belonging to a triplet state, 

£t — £a + e& + Cafe — ^ab, 

and one eigenvalue belonging to a singlet state, 

e s = e a + e + Cab + J a b ■ 



(5) 



(6) 



(7) 



(8) 



(9) 



(10) 



The exchange integral can be shown to be positive: if we take /(r) = a (r)0 b (r) 
and perform Fourier transforms, we have 

o2 



J ah = I dY l f*(r 1 ) [dr 2 - — r/(r 2 ) 

J j \Ti — r 2 



(2tt) 3 / 2 
47re 2 



rfke 



/ — ik'ri f* (\l 



/*(k' 



(2tt) 3 
rfk|/(k)| 2 



47re 2 



A: 2 



> . 



(11) 



Because J a b is greater than zero, the triplet is lower in energy than the singlet; 
i.e. in the lowest state the spins are parallel, which corresponds to a ferromagnetic 
situation. 
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The action of the Hamiltonian can be expressed by spin operators. We have 

^ ^ , ,„ „ i 1 for the triplet , , 

2SiS 2 + l/2 = (Si + S 2 ) 2 -l = ^ P . (12) 

[ — 1 for the singlet 

The action of the triplet and singlet can then be expressed by a single Hamiltonian 

U = - ^ y ^(2S 1 S 2 + l -) = const - 2J a6 S 1 S 2 . (13) 
Generalizing the exchange interaction to a lattice, one may write 

W = "^E-Wi . (14) 

This is the most familiar form of the Heisenberg exchange, where % and j represent 
lattice site indices and the factor \ is introduced by convention. 

2.2. Direct exchange with non-orthogonal states (antiferromagnetism) 

Antiparallel spin alignment (antiferromagnetism) occurs in a two-center system like 
a hydrogen molecule in the Heitler-London approximation. Consider two hydrogen 
atoms centred at R a and Rf, respectively, with a Hamiltonian in which each elec- 
tron feels both protons and the electron-electron and proton-proton interactions are 
included 

7~t = H atom (r 1 — R ) + H atom (r 2 — R&) 

e 2 e 2 e 2 e 2 

—. r- 1 - + - - + - r . (15) 

Iri-Rftl |r 2 -R a | |r x — r 2 | |R a - «*| 
Each electron occupies a separate ls-orbital centred on one of the atoms. In the 
simplest approximation, the low-lying states are assumed to be described by four 
configurations with spins T1\tl>l1\ll- We denote the orbital basis functions at 
atoms a and b as </> a (r) and </>&(r). They are in general non-orthogonal: 

J 0I(r)0 6 (r)rfr = 1^0. (16) 

Each orbital function is associated with one of two spin functions, xt or Xi- The 
Hamiltonian is diagonalized in the subspace of the following four normalized spin- 
coupled functions corresponding to triplet and singlet states tp(S, S z ): 

VKM) = , = XT( g i)XT( s 2)[0a(ri)0b(r 2 ) - 0„(r 2 )0 6 (ri)] , 

= / = o fXT( g i)Xl( s 2) +X|(gi)XT( g 2)][0a(ri)0 fc (r 2 ) -0 a (r 2 )0 b (r!)] , 
/2(1 - P) 



^(1,-1) = i 1 X|(gi)X|(g2)[0a(ri)0 b (r 2 ) - a (r 2 )0 fc (ri)] , 
^2(1-/2) 

^(°'°) = T, A x\{s x )xi{s 2 ) - X|(gi)XT( g 2)][0a(ri)0 fc (r 2 ) + a (r 2 )0 6 (ri)] 

2^(1 + / 2 ) 



(17) 



7 



The triplet energy may again be written in terms of Coulomb and exchange integrals: 

e t = 1)) = 2e atom + Ga [~ I { ; h . (18) 

Here, 2e a t om comes from the one-electron part of the Hamiltonian. The Coulomb 
integral (containing terms where the electron belonging to one nucleus feels the 
attraction of the other nucleus) is 




and the exchange integral is 

I ab = [ dr 1 [ rfr 2 0*(ri)0 f) (ri) 1 — r <^(r 2 )0 a (r 2 ) 

J J \ti — r 2 | 

2 2 

-/ / rfri 6 M tJMti) ~ I I dv 2 . \ <^(r 2 )0 a (r 2 ) . (20) 

J \Ti — ttb\ J 1 1*2 — Ho | 

The singlet eigenenergy is 

e s = (V»(0,0)|W|V(0,0)> = 2e atom + ° ab ± [ ab . (21) 

1 + r 

The singlet-triplet splitting is 

e t -e s = 2 l2C i ab _- 4 hb . (22) 
As in the previous subsection, an effective Hamiltonian may be defined as 

TC = const + J 12 SiS 2 , (23) 

with 

Ji 2 = 2 ^^ 6 _~/ a6 . (24) 

Without the Coulomb term, one has ferromagnetic coupling; with a sufficiently large 
overlap, the effective exchange coupling becomes antiferromagnetic ( J 12 > 0); the 
ground-state of the hydrogen molecule is a singlet state. Generalizing to a many- 
electron system, one has the Heisenberg model for an antiferromagnetic lattice. 

Including ionic configurations where both electrons can sit on one or the other 
atom leads to a hopping mechanism for the electrons (kinetic exchange), which 
supports antiferromagnetic coupling, see e.g.(jS], p. 60). 

An antiferromagnetic Heisenberg model is also obtained from the Hubbard model 
in the strong coupling limit, or from indirect exchange mechanisms like the RKKY 
scheme, leading to an effective Heisenberg model in second-order perturbation the- 
ory. Also, super-exchange or double exchange lead to Heisenberg like terms or even 
to biquadratic terms, see e.g. 
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In the present article, we do not try to give a better justification of the Heisenberg 
model. Rather, we consider it a phenomenological model that proves to be successful 
in describing many experimental data when its parameters are fitted. A Heisenberg 
model is adequate when the spins are localized (e.g. in the rare earth elements). It 
should also be applicable to 3d-transition metal band magnets because the magnetic 
moments are quasi-localized when integrating over microscopically calculated spin 
densities. One also sees in experiments on bulk transition-metal ferromagnets that 
the magnetization follows the Bloch T 3 / 2 law at low temperatures. Above the Curie 
temperature, one observes a Curie- Weiss behaviour of the magnetic susceptibility. 
Both features follow from a Heisenberg type model. 
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3. Basic equations of the Green's function formalism 

In this section we place together the essential definitions and derivations of the 
Green's function formalism which are necessary to understand the following article 
without frequent recourse to the literature. For further details of the basic features 
of the Green's function formalism as it is used in the present review, we recommend 
the article pQ and the books |3] and [Sj. 

The double-time Green's functions (GF's), as they are exclusively used in the 
present article, are defined in Section 3.1 and their equations of motion are given 
in Section 3.2. In Section 3.3 we discuss the eigenvector method for determining 
the GF's. Once the GF's are known, the corresponding correlation functions (ther- 
modynamic expectation values) are determined by the standard spectral theorem, 
where, in general, commutator and, in the case of zero eigenvalues of the equation- 
of-motion matrix, anti-commutator GF's have to be used. A proof of the standard 
spectral theorem is given in Section 3.4. In Section 3.5 we discuss the singular value 
decomposition of the equation-of-motion matrix and show how a transformation can 
be found to eliminate the null-space, obviating the need for the anti-commutator 
GF. This procedure is necessary whenever the quantities associated with the null- 
space are momentum-dependent because the standard spectral theorem fails in this 
case. In Section 3.6 we show that there is no advantage in starting the calculations 
with the anti-commutator GF instead of the commutator GF. In Section 3.7 we show 
how the intrinsic energy, the specific heat and the free energy can be calculated with 
Green's function theory (GFT). 

3.1. Definition of the double-time Green's function 

Because we deal later with multi-dimensional problems, we prefer to work with a 
vector of Green's functions having components characterised by the index a: 

Gf j>v (t-t') = {{Amm')))n = -m-mA^B^x) . (25) 

Throughout the paper we deal exclusively with such double-time GF's. In principle, 
either the commutator (77 = —1) or anticommutator (77 = +1) of the Heisenberg 
operators Ai(t) a and Bj(t') can be used (but see Section 3.6); i and j are lattice site 
indices. 

The operators obey the Heisenberg equation of motion, e.g. 

A2(t) = e iHt Afe- iHt , A? = -i[Af, K\_ x . (26) 
Here H is the Hamiltonian under consideration. We set ft = 1 throughout the paper. 
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For magnetic films, the Af are spin operators obeying the usual commutator 
rules. One has, for instance (see Section 4.2.1), 

Af = {St, ST, St), (27) 
Bj = (S*) m (Sr) n with m + n < 2S + 1 (m > 0, n > 1, integer). 

In eqn f!25|) . i is the imaginary unit (when it is not an index) and the step function 
Q(t — f) is defined as 

f 1 for t > t' 

0f-f =i (28) 
[0 for £ < f . 

The double brackets ((Ai(t) a ; Bj(t'))) n are an alternative notation for the Green's 
functions Gfj Jt—t'). Single brackets denote correlation functions, e.g. ([Af (t), Bjit')]^), 
which are thermodynamic expectation values 

(•••) = \ £(n|e-^..» = ^Tr(e"^...) , (29) 

where 

Z = J2(n\e~ pH \n) = Tr(e- f3H ) (30) 

n 

is the partition function with (3 = 1/(A;^T), T the temperature and ks the Boltz- 
mann constant. 

Usually, it is more convenient to work with the Fourier transforms of the Green's 
functions in energy space, 

r°° 

Gfj(u) = / dit-t'^t-t'y^, 

' — oo 

i- 00 du 

'-oo 2-7T 

and in momentum space, 

ij 

G%{u>) = ^E^(^)e- ik(R - Rj) , 
iV k 



GUt-t') = f ^Gf^)e~^'\ (31) 



with &: = iVe^-^l 



k 



and <W = i^e^. (32) 

i 

Here the Rj are the lattice site positions and N is the number of lattice sites. 
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3.2. The equations of motion 

The Green's function vector has to be determined by its equation of motion. This 
is obtained by taking the time derivative of equation (|23jl . 

l Ft G? ^ {t ~ f) = 5{t ~ WW®' B ^ + «^ ' B M))v > (33) 

where the Heisenberg equation (}2l)j) has been used together with j^Q(t—t') = 8(t—t'). 
Eqn (|33J) is a differential equation for determining the Green's functions. Because it 
is more convenient to work with algebraic equations, one usually performs a Fourier 
transform to energy space (j31|) . characterized by the index u: 

= ([A?, Bj] v ) + (([A?, #]_; S,-))^ . (34) 

Observe that on the right-hand side a higher-order Green's function arises which 
leads to another equation of motion having even higher-order Green's functions and 
so on. In this way, an exact infinite hierarchy of equations of motion is generated. 
Only in rare cases does this hierarchy terminate automatically. Usually, one has to 
terminate the hierarchy somewhere in order to obtain a solvable closed system of 
equations: the Green's function of some specified order must be factored in such a 
way as to contain only Green's functions which already exist in the hierarchy up to 
the cut-off. This factorization is called the decoupling procedure and is the essential 
and most severe approximation in GF theory. Except for a few cases, it can be 
justified only by its success. 

Very often one works with the lowest-order equation (J34)) only. The decoupling 
consists in this case in factoring the GF: 

(<[A^]-;^>>^££rf {{a^b,))^ . (35) 

l (3 

The right-hand side now has only GF's which are of the same order as those already 
present. In this way one arrives at a closed system of equations of motion. The 
matrix T°f is in general unsymmetric. 

Inclusion of the second-order equation of motion would require a decoupling of 
the double-commutator GF (([[Af , H\, H]; Bj)), etc. 

For periodic lattice structures, the equations of motion are simplified by a Fourier 
transformation to momentum space (}3*2*j) . which eliminates the lattice site indices. 
The equations of motion in compact matrix notation are then 

(col - r)G„ = A, , (36) 

where is the inhomogeneity vector with components A® = ([A a ,B] v ), and 1 is 
the unit matrix. 
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From Kramers rule one sees that the GF's have a pole structure with the eigen- 
values of the matrix F as poles. In many applications, use is made only of these 
eigenvalues but we show in the next subsection that it is of great advantage to 
use the eigenvectors of this matrix as well, especially in treating multi-dimensional 
problems. 

3.3. The eigenvector method and the standard spectral theorem 

In this section, we show how to take advantage of the eigenvectors of the matrix 
r in transforming the GF's to a new set of GF's each having but a single pole. 
This is particularly important in treating degenerate eigenvalues of T because each 
eigenvalue can be associated with a definite (transformed) GF. Also, the extra cost 
of finding the eigenvectors is more than compensated by avoiding the effort of cal- 
culating determinants in a Kramers-like treatment and by the clarity gained in the 
formulation. We shall use the notation of reference 0. 
The first step is to diagonalize the matrix T 

LTR = n, (37) 

where f2 is the diagonal matrix of N eigenvalues, u T (r = 1, N), No of which are 
zero and (N — N ) are non-zero. The occurence of zero eigenvalues is not a rare 
case: they arise as a consequence of the spin algebra of certain of the GF's and are 
to be expected. The matrix R contains the right eigenvectors as columns and its 
inverse L = R 1 contains the left eigenvectors as rows. L is constructed such that 
LR = 1. We assume that the eigenvectors span the whole space so that it is also 
true that RL = 1. 

We now define new vectors by multiplying the original vectors with L: 

g v = LG V and A v = LA V . (38) 

Multiplying equation (|36|) from the left by L and inserting 1 = RL leads to 

(wl - fl)g v = Ar,. (39) 

From this equation we see at once that each of the components r of this Green's 
function vector has but a single pole (!) 

{S„)r = ^ • (40) 

This allows a direct application of the standard spectral theorem (see e.g. 
for its proof, see Section 3.4) to each component of the Greens's function vector 
separately. The spectral theorem relates the correlation vector 

C = LC k = L(BA) k (41) 
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to the Green's function vector. The index k indicates that we work in momentum 
space. Explicitly, 



Cr = ± Mm r ^(g.t" = (A^_ (42) 

where eqn (HUJ) and —r^ = £ =F «7r5(a;) have been used. 

In general, we use the commutator GF's, in which case the inhomogeneities 
A n =-i are independent of the momentum k, whereas the A v=+ i are not. Using the 
anti-commutator GF's (rj = +1) leads to problems connected with this k-dependence 
(see Section 3.6). The commutator GF's, on the other hand, lead to problems with 
zero eigenvalues of the equation of motion matrix T because there are then zeroes in 
the denominator of eqn (|42p. In this case, the correlation vector must be split into 
two components C\ and C° o belonging to non-zero and zero eigenvalues respectively. 
We then have 

^ = £=r < 43 > 

where u T ^ 0. 

For the correlation vector belonging to zero eigenvalues, the anti-commutator 
GF is required (for a proof, see Section 3.4): 

(C°) ro = lim -a;(^ =+1 ) T0 
oj— >o z 



2 u - (uo TO = 0) 2 V + ; ° 
= l(L (A_ 1 + 2C k )k = (L C k ). • (44) 

Here, the relation A + i = A_i + 2Ck has been used together with the fact that the 
commutator GF is regular at the origin (see eqn (J7DJ)). By multiplying eqn (J36j) with 
L°, using L°r = and taking the limit u — > one obtains 

limL°(cul-r)G_i = L°A_ 1 = 0. (45) 

We call this the regularity condition. 

We now partition all quantities with respect to the non-zero and zero eigenvalue 
space 

, i nx I Ll \ I Cl = PV^-l \ , , 

R = R^ L = ,C= , 46 

1 ' V L ° / \ C° = L°C k ) K J 

where S 1 is a diagonal (N — N ) x (N — N ) matrix with elements l/(e f3wr — 1) on 
the diagonal (uj t ^ 0). 

The original correlation vector in momentum space is then 

RC = R'R ) I ^ () | \i : I. \ , +R"L ,, C k . ( 171 



c 
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We are interested in the diagonal correlations C (without the index k) in configu- 
ration space 

C = h E C k = / ^C k , (48) 

k 

where the integration is over the first Brillouin zone. This leads to a set of integral 
equations for the components Cj (i=l,...,N) which have to be solved self-consistently. 
If the factor R°L° is momentum independent, one can take it outside the integration 
in the second term of eqn (}4T|) to get the Cj components explicitly: 

. N-N N N N 

Ci = / dk( E E /«'/;,/-;/(•! + E E 4 L ?A m 

3=1 1=1 j=ll=l 

The components Cj are obtained by iterating on the C» until Eqn ()49j) is satisfied. 

If R°L° is momentum-dependent, the standard procedure fails because one can- 
not take R°L° outside the integration. Instead, one needs a more complicated 
procedure that relies on the singular value decomposition of the T-matrix (see Sec- 
tion 3.5). This leads to a formulation of the spectral theorem in which the null-space 
is eliminated and only the commutator GF is needed, obviating the use of the anti- 
commutator GF. 
oo 

Equation (|47[) can also be derived without the anticommutator GF in the following simple 

way: 

Start with the spectral theorem for the commutator GF (|42|) with rj = — 1 

C k = R£'LA_ 1 , (50) 
and make use of the decomposition (|46[) and £q = l/(e^° — 1) = oo. Then 

C k = (R 1 R°) ( ^ ° Q j ( ^ j A_x = Rif^A-i + K°£°L°A^. (51) 

The second term is undetermined because £ °L°A_i has the indeterminate form oo x 0, 
see ean (|45|) . We get around this by multiplying the last equation from the left by R°L°: 

R°L°C U = R^R^l/A-i + R°L R £ L°A_i = R°£ Q L°A_i. (52) 

The last term is obtained because L°R 1 = and R L°R°L = R°L°. Thus the term 
R°£ L°A _i in eqn (|5T|) can be replaced by R°L C k , which completes the proof of eqn 

gZJ). 

oo 

3.4. The proof of the standard spectral theorem 

The spectral theorem is the relation of greatest importance for the Green's func- 
tion formalism because it allows the calculation of the desired observables (or more 
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generally the correlation functions) from the corresponding Green's functions. Al- 
though its proof can be found in text books (e.g. Ej), we reproduce it here for 
the convenience of the reader. 

Considering one component of the GF vector (|25|) (for brevity we leave out the 
index a), we introduce the spectral function Sij jV (t — t') by 

G ij>n (t - f ) = -iQ(t - t')2irS ijiV (t - 0, (53) 

where, by comparing with eqn ()25|). 

S lhV (t-t') = ^([Mt^B^X) = ^-{AiifyBjit') +r}B j (t')A i (t)) . (54) 

Inserting a complete set of eigenstates (H\m) = uj m \m)) yields the following spectral 
representations for the correlations: 



{A t (t)B,(lf)) = - V{n|B^|m){m| J 4^|I^)e-' 3 ""e s( ""-""■ ) e- i( "•-""' , "- , ' , , (55) 



nm 

1 





\m)(m 


\A\ 


\B 3 


\m)(m 





(Bj{t')Ai{t)) = ^^(n|5 j |m)(m|A»e-^"e- 4{ ^-" m){t - t ' ) , (56) 

Zj 



nm 



and the spectral function, 

Sij, v (t - f) = — ^(n| J B i |m)(m|A»e- /3 ^(e /3( "' l - Wm) + ^e"^'— )(*-*'), (57) 

whose Fourier transform to energy space is 

S ijiV (u) = jY,( n \ B M)HMn)e-^(e^ + r ] )5(iu-(uj n -u Jm )). (58) 

^ nm 

A relation between the energy representations of Sij tV (uj) and Gij tV (uj) is derived by 
inserting in 

G ij:V (uj) = -2ni / d(t - t'y^*-* >0(t - t')S lhV {t - t') (59) 



the following representation for the step function 

Q(t-t') = - dx 6 -——, (60) 

Z7T J-oo X + IT! 



and the Fourier transform of Sij, n (t — t'): 



/oo roo 1 1 poo 

duj' / dx— — / d(t - t')e<^ 
-oo J-oo X + lT]2nJ-oo 



" d J S J2fk^l . (61) 

oo UJ — UJ + 11] 



With 



1 1 

Gij.^uj + i8) - G ijiV {uj -iS)= I duj' S ijyV (uj')( — — -) (62) 

3 UJ — uj' + 10 UJ — UJ — 10 
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and 

1 T—t = P—^ T i*6(u - uj') (63) 

uj — uj' ± id uj — uj' 

it follows that 

S ijiV (oj) = \im-^-(G ijri (u + i5) - G ijtV (uj - i5)) . (64) 

We can also see that 

(BjWMt)) = -^—S^e-^-^ (65) 

by inserting equation ()58|) in this equation and comparing with ()56|) . 
Together with equation this yields 

(b j( o^w> = J™ f /" ^" + ^- g ^"- tf ) e -^) . (66) 



This is nothing else then equation ()42|) with t = t' after a Fourier transformation to 
momentum space. It is valid for rj = ±1. 

For 77 = — 1, this expression diverges in the limit uj — > and it is necessary to 
use eqn (jUJ). This was first pointed out in reference fTUj . see also [TT], and can be 
seen by decomposing the spectral function (J58|) into two terms referring to uj n 7^ c<j m 
and u n = uj m respectively 

S ij>v (u) = S ij>v U^ m + (1 + rj)C%8{u). (67) 

Inserting this in eqn (|61|) and taking the limit uj — > of ujGij tV (uj) one finds 

lim a;Gj, „(u;) = lim / cL/ f — — — h - — ^LUl — l J-S — \\ 

= + (l+77)C°. . (68) 

From this expression, we see that the quantity Cfj is determined by the anti- 
commutator GF (rj = +1) 

G% = \\\^ujG im= + x {yj) , (69) 

whose Fourier transform to momentum space is equation ()44j) . This completes the 
proof of the standard spectral theorem. 

From equation ()68|) an important analytical property follows: the commutator 
Green's function (rj = —1) is regular at the origin, 

lim^G ijfl= _i = 0, (70) 

a fact which is necessary to derive the regularity condition (|43|) . The anti-commutator 
Green's function has a first order pole at uj — 0. 
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3.5. The singular value decomposition of T and its consequences 

In this section we show that the singular value decomposition of the equation-of- 
motion matrix T obviates the need to use the anti-commutator GF when zero eigen- 
values occur; the commutator GF suffices. 

The standard spectral theorem is of practical use only if the quantity R°L° in eqn 
(147)1 is momentum independent, because only then can one arrive at an equation that 
can be solved by iteration (see ([49)1 . If R°L° depends on momentum, the standard 
procedure fails because equation (|4Tjl is of the form 

(l-R°L°)C k = R 1 ^ 1 L 1 A_ 1 . (71) 

The term (1 — R°L°) is idempotent and therefore has no inverse; hence, one cannot 
solve for Ck- This arises for instance for the reorientation of the magnetization using 
exchange anisotropies, see Section 4.2.3. 

oo 

An idempotent operator P has no inverse. 

Proof: assume the existence of an inverse: P~ 1 P = 1 and idempotence P = P 2 , 
then P _1 P 2 = 1, implying P = 1, which is a contradiction. 

oo 

The singular value decomposition (SVD) offers a way out of this situation by 
providing a transformation that eliminates the null-space; in effect, it defines a 
smaller number of Green's functions whose associated equation of motion matrix, 
7, has no zero eigenvalues, thus dispensing with the anti-commutator GF as well as 
reducing the number of equations. 

The singular value decomposition states that ... "any M x N matrix A whose 
number of rows M is greater or equal to its number of columns, can be written as 
the product of an M x N column-orthogonal matrix U, an N x iV diagonal matrix 
W with positive or zero elements and the transpose of an N x iV orthogonal matrix 

v..., m 

The equation-of-motion matrix can therefore be decomposed as 

r = UWV = uwv. (72) 

where U and V are orthogonal matrices (UU = 1, VV = 1) and W is a diagonal 
matrix with singular values on the diagonal. U, V and W can be determined very 
efficiently numerically [^j. The matrices U and V can also be obtained by diago- 
nalising IT or IT respectively. The singular values are the positive square roots of 
the eigenvalues of these matrices: 

vfrv = VVWUUWVV = W 2 , (73) 

urfu = uuwvvwuu = w 2 . 



18 



If r has zero eigenvalues, it has the same number of zero singular values. The matrix 
r is also given by uwv, where u and v are obtained from U and V by omitting 
columns corresponding to singular values zero, u and v are again orthogonal matri- 
ces (uu = 1, vv = 1). Note that vv is a projector onto the non-null-space and v v 
a projector onto the null-space (vv + v v = 1). The matrix w is diagonal having 
positive singular values on the diagonal. 

To eliminate the null-space, it suffices to use the following transformations: 

7 = vTv, 
g = vG, 
a = vA, 

c = vC k . (74) 

Multiplying eqn (J36|) by vv = 1 and T = uw(vv)v = (uwv)vv = Tvv one obtains 

v(uj1 — Tvv)G = vA_i, 

(ul — vTv)vG = vA_i, 

(ul - 7 )g = a. 

Now we diagonalize 7 

l7r = lj 1 , where 1 = L 1 v and r = vR 1 . (75) 

7 is a reduced matrix with the same non-zero eigenvalues u 1 as the original matrix 
r. Since there are now no zero eigenvalues, we can apply the spectral theorem with 
respect to the non-null-space: 

c = rSHa, (76) 

where S l is the matrix occurring in eqn ([47)1 . A Fourier transformation to configu- 
ration space yields the self-consistency equations (analogous to eqn 



= J dk(rf 1 lvA_ 1 - vC k ). (77) 

Again, this can be solved for the correlations in configuration space C if one can 
find a row- vector Vj which is k-independent, i. e. 

J dk Vi C k = Vj J dk C k = VjC. (78) 

This equation may be supplemented by the regularity condition 

lim u (c<j1 — uwv)G = uoA_! = 0. (79) 

This is because uqu = and because the commutator GF is regular at the origin. 
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One may be tempted to object that eqn (J78j) is no improvement over eqn (J49j) 
because, in both cases, it is the k dependence of a term containing Ck that creates a 
problem. In practice, however, it is much better to use SVD because diagonalization 
of the full matrix T to get R° and L° is fraught with numerical difficulties when 
there are non-zero eigenvalues which are very small. Furthermore, the vectors R° 
and L° are non-orthogonal, whereas the projector vv onto the non-null-space is built 
from orthogonal vectors - this makes it easier in practice to find a row vector Vi that 
is independent of the momentum k. This search is technically complicated, and for 
a more detailed description, we refer the reader to Ref. [T3|. Here, we give a recipe 
for seeking for appropriate Vj. 

The row vectors in v are determined numerically and are unique up to a sign 
change or, for degenerate singular values, up to an orthogonal transformation of 
the degenerate vectors. In order to distinguish among the row vectors of v, it is 
very helpful if they are suitably labelled; e.g. they can often be characterized by a 
layer index or a sublattice index. The following procedure is useful: decompose the 
T-matrix into a reference matrix and the rest, 

r = T re f + T rest , (80) 

where T re f has a block structure determined by the chosen labels. With a singular 
value decomposition, 

Tre/ = U r e/W r e/V re j, (81) 

one can define a block-label operator 

N B 

p ° P : =E V -/(^(*)V re/ (*), (82) 
i=i 

with L(i) = Nb — i + In the basis of the singular vectors v (and analogously for 
vo), we define a matrix 

N B 

P = vP op v = J2 vV re/ (i) V^W VM0V re/ (z)v = SS (83) 
i=i 

with 

S = [y/Z(tjVref(l) © - © L(N B )V ref (N B )]v , (84) 

where © defines the direct sum. 

Now the singular value decomposition of S, 

S = LyZ, (85) 

furnishes a matrix Z that diagonalizes SS: 

SS = ZyLLyZ = Zy 2 Z, (86) 
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where y 2 L(i), which labels the blocks. To each block-label belongs a labelled 
vector 

v L = Zv, (L = labelling) (87) 

which is the desired result. 

A further difficulty is connected with the fact that the computed v will not 
necessarily be continuous, even if the elements of the T-matrix are changed contin- 
uously (e.g. by varying the momentum k on which they depend); i.e. vectors at 
neighbouring values of k can have arbitrary phases. This difficulty is overcome by 
a smoothing procedure, which consists of the following steps: 

(1) Create well-behaved reference vectors V re j(r = l,..,N r ) in the momentum 
range of the first Brillouin zone for the vectors V° at k and V 1 at k±, etc. by 
overlaps as in the labelling procedure. 

(2) Interpolate the reference vectors at each k 

V ref (k) = w{V ref (ki) + w h ~V ref (k h ) (88) 

with 

o /7T . k — ki . \ . , 

wi = cos {^ kh - ki ' w h = l-wi. (89) 

(3) Orthonormalize the reference vectors 

y = V re /V re /, 

a = fyT, 

y-1/2 = TX -1/2 T ^ 
V re/ = y- 1/2 V ref - (90) 

We now have reference vectors for the non-null and the null-space: V re j = (v re j, v re j). 

(4) Match the untreated (or, if necessary, labelled) vectors v to the orthonormal- 
ized reference vectors vy e /. This is done by seeking a transformation Q that rotates 
the target (original) vectors among themselves to achieve the best match 

vs = Qv. (S = smoothed) (91) 

Q is found by a SVD of the overlap matrix of the reference vectors with the 
target vectors 

S = v re/ v = £xl. (92) 

Here x is a diagonal matrix of the singular values of the overlap matrix which are 
close to 1 by construction (x ~ 1). The desired transformation matrix is 

Q = CZ, (93) 
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which is a rotation matrix because 



QQ = ZtLZ = ZZ = \. (94) 

The overlap matrix of the reference vectors with the smoothed vectors is close to 
the unit matrix because the phases of the new vectors have been fixed by Q: 

Vre/Vs = v re /vQ = SZC = CxZZC = £x£ ~ 1. (95) 

To summarize, the untreated vectors v of the original problem can be labelled and 
smoothed by the tranformation 

vls = Qv L = £ZZv. (96) 

In practice, some of the row vectors of this transformation matrix turn out to be 
momentum-independent and can be used in solving equation (J77J). 

The procedure described above was successfully applied to Heisenberg multi- 
layers with exchange anisotropies, see Section 4.2.3 and to coupled ferro- and an- 
tiferromagnetic layers, see Section 4.3.2. We stress once more that the standard 
spectral theorem fails in these cases. 

3.6. No advantage to using the anti-commutator instead of the commu- 
tator Green's function 

We begin with the simplest case of a Green's function G v which has but a single 
pole and an inhomogeneity A v : 

G v = ((A; B)), 

A v = ([A, B] v ), (97) 

i.e. 

The corresponding correlations in momentum and configuration space are 

C k = (BA), (99) 

c = -yc k . 

Applying the spectral theorem gives 

C k = s^— ■ (100) 
+ r] v ' 

Note that A + i(k) = A_i + 2Ck and A + i(k) depends on k, whereas A_i does not. 
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The commutator (77 = —1) GF yields the correlation in configuration space 

< 101 > 

whereas the anti-commutator (77 = +1) GF leads to 

civMav^^, (102) 

N ^ e ^ + 1 N ^ + 1 v ; 

which cannot be solved because C k is unknown. Putting eqn (jl(J(Jj) with 77 = — 1 
into eqn ()102|) leads again to equation (jlOlj) . 

c- 1 y A - 1 + 2g fe - 1 v A - 1 rin^ 

AT V e/*-* + 1 AT V - 1 ' 1 ' 

k k 

which can be solved self-consistently. This shows that there is no advantage in 
starting the calculation with the anti-commutator GF. 

One can show this more generally with the eigenvector method of Section 3.3, 
see |2Zj: starting with the anti-commutator formulation, the spectral theorem yields 

C k = R£LA rj=+1 , (104) 

where £ is a diagonal matrix with elements 5^ = 8ij(e^ U3i + l) -1 and A v=+ i depends 
on the momentum k, preventing a direct use of this equation. 
Because A )?=+1 =A_ 1 + 2C k , 

C k = R£L(A_x + 2C k ) (105) 

C k = (1 - 2R£L)~ 1 R£LA_ 1 . (106) 
(R£L) _1 = L-^R- 1 = R£ X L (107) 



or 



Introducing 



in <nm . 

C k = (R(l - 28 )L) _1 R£LA_i = R(l - 2£)~ 1 £LA_ 1 = R£LA_ 1; (108) 

where = 5ijS„/(l — 2Su) = 5 i j(e^ i — l) -1 . This is still of no use because of the 
zero eigenvalues. But we have shown in Section 3.3 that the term R°L° remedies 
this: 

C k = R 1 £ 1 L 1 A_ 1 + R°L°C k (109) 
which is equation (J47|) . where S 1 is the matrix £ leaving out the diverging terms. 
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3.7. The intrinsic energy, the specific heat and the free energy 

The intrinsic energy is the thermodynamic expectation value of the underlying 
Hamiltonian 

E={H) = NE i} (110) 

where E\ is the intrinsic energy per lattice site and N is the number of lattice sites. 
The specific heat at constant volume is obtained by differentiating the intrinsic 
energy with respect to the temperature 

dE „ dE 

cv = If = Tp ■ < m > 

The free energy is obtained by integrating over the intrinsic energy 

F(T) = E(0) -T f T dT' E ( T ^ ~ E ^ . (112) 
Jo T' 

oo 

Proof of this formula: 

From F = E — TS and S = — 4y | y one has 

, d F(T) 

E(T) = —T 113) 

dT T K J 

from which one obtains eqn (|112|) by integration. Differentiating (|112|) gives 

oo 

In order to see how the intrinsic energy per lattice site can be calculated explicitly, 
consider the quantity 

B? A = (A i [C i ,H]_), (114) 

where Ai and Cj are the spin operators necessary for constructing equations of 
motion for those Green's functions from which the moments ((S z ) n ) (n = 1, ...,2S) 
are calculated. 5 is the spin quantum number. The quantity ()114|) can on one hand 
be related to the relevant Green's functions and on the other hand be calculated 
explicitly by evaluating the commutator. This leads to a set of equations from which, 
together with equation (|11U|) . the intrinsic energy can be calculated. 

The connection to the Green's function results from the spectral theorem: 

d_ 
dt 

limiV-f _^_( G ^( W + iS) - G°> A (u - tf)) e -M*-0| 
s^o N V2W eft*- 1 



Bf> A = (A l {C i ,H}„)= l -(A i (t')C l (t))\t=t> 
. d 1 % r duj 



dt s^o N ^ 2tt J eP"- V k v ; k v 1 

k 

= ft^s/ 3^l( G ^ <<" + "> " ^ ~ "» ' < 115 > 

k 

In Appendix 7.1, we treat explicitly the cases for spin 5 = 1/2 and 5 = 1 for 
a Heisenberg Hamiltonian with an external field and a single-ion anisotropy. For 
S=l/2, one needs Ai = S~ and Cj = 5 + ; for S=l, one needs a) Ai = S~ and 
d = 5+ and b) Ai = 5^ and C t = (25? - 1)5+. 
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4. The GF formalism for Heisenberg films 

This chapter starts in Section 4.1 with the example of a spin S = 1/2 ferromagnetic 
Heisenberg monolayer in a magnetic field This is an exercise in applying the 
GF formalism in a simple case. The Tyablicov (RPA) and Callen decouplings are 
introduced, the limit of mean field theory (MFT) is discussed, the Mermin- Wagner 
theorem is proved for this case, and the effective (temperature-dependent) single-ion 
anisotropy is calculated by thermodynamic perturbation theory. 

In Section 4.2, ferromagnetic Heisenberg films with anisotropics and general spin 
S are treated. For the single-ion anisotropy, the Anderson- Callen decoupling is 
used. The exchange anisotropy is treated by a generalized Tyablikov decoupling. 
Susceptibilities are calculated. It is also shown how the single-ion anisotropy can 
be treated exactly. As a further application, it is shown that spin waves are very 
important for treating a trilayer in which two ferromagnets are separated by a non- 
magnetic layer. Finally, the temperature dependence of the interlayer coupling is 
discussed. 

Section 4.3 deals with a unified treatment of ferromagnetic (FM), antiferromag- 
netic (AFM) and coupled ferromagnetic-antiferromagnetic (FM-AFM) Heisenberg 
films. 

4.1. The ferromagnetic Heisenberg monolayer in a magnetic field 

We choose this example because it illustrates the GF formalism in a simple case 
and allows the validity of the different approximations within the formalism to be 
checked against 'exact' Quantum Monte Carlo (QMC) calculations |15j . 

The Heisenberg Hamiltonian for a ferromagnetic monolayer in a magnetic field 

is 



Here Jki is the exchange interaction strength, k and / are lattice site indices, and 
< kl > means summation over nearest neighbours only. The magnetic field B is 
assumed to be in the z-direction perpendicular to the film xy-plane . The second 
line of eqn (fTTBj) is obtained with the usual definition S k = 5f ± iS% in terms of the 
components of the spin operators. 

For spin S — 1/2, the magnetization is obtained from the relation 



H 



^/SfcS/ — b y sf 



<kl> l 





(Sf) = 1/2 - (SrSf) 
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and the correlation (S { S^~) is determined via the spectral theorem from the com- 
mutator GF 

G^ 1 (u) = ((St;Sr)). (118) 
The GF is determined from the equation of motion in energy space 

co((S+;Sr)) = ([S+, Sr)) + {([S+, H].;Sr)) . (119) 

Using spin commutator relations, one obtains 

[S+, H}_ = BSf - £ MSfSf - SfSf). (120) 

i 

The equation of motion is then 

(u - B)((St; Sr)) = 2<Sf)<% - £ J«(«5?#; Sj)) - ((SfSf] Sj))) , (121) 

i 

which is exact as it stands but, in order to use the equation, the higher-order Green's 
functions on the right hand side must be decoupled. 

4.1.1. The Tyablikov (RPA)-decoupling 

This decoupling, introduced by Tyablikov [IB] , is often called the random phase 
approximation (RPA) because it is equivalent to that approximation in other areas 
of physics. It consists in factoring the higher-order Green's functions: 

* (Sfi((S?;Sj-)) = (SfiG Ij , 
((S?St;Sr)) ~ (Sn((S+;S.)) = (SnG ir (122) 

There is no a priori justification for this factorization but it has turned out to be 
successful, also in other areas of physics where the resulting equations can be derived 
with methods different from Green's function theory. In the present context, the 
quality of this approximation can be checked against 'exact' QMC results [T5] . see 
Section 4.1.5. 

For a ferromagnet, there is translational invariance for the magnetization at 
different lattice sites: (Sf) = (Sf) = (S z ). After the decoupling, the equation of 
motion is 

{u-B- (S z ) Ju)GiM + (S z ) E JuGiM = 2{S z )5 ir (123) 
i i 

A Fourier transform to momentum space (J32)) yields 

(u-B- (S Z )(J - J k ))G k (u) = 2(S Z ), (124) 
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and the Green's has the pole structure 

Gk(«) = 2{S 7pa > (125) 

with the dispersion relation 

u* PA = B + (S*)(J - ^k). (126) 

For a square lattice with the number of nearest neighbours 2: = 4 and a lattice 
constant unity, one has 

Jo = ^E^ i(k= ° )(Rl ^ ) = ^ = 4J, 

ij 

J k = ^E^ e4k(R '" Rj) = 2J(cOSfc e + COS fcj,). (127) 

ij 

Applying the spectral theorem (|42|) - there is no zero eigenvalue - and performing 
the ^-integration with the relation 

1 „ 1 



=F m6(uj - wt) (128) 



yields for the magnetization of the spin 5=1/2 monolayer 

1 1 i roo Gf k (u; + i5) - G k (^ - id) 



lim — — f 

s^o N , 2n J- 



2 tf-o iV Y 2vr 7-oo e#" - 1 

1 1 y, 2(S Z ) 



2 N Y e^k - 1 

= l~~ 2 Tdk x Tdk y -ff (129) 
2 7r jo jo k — 1 

where the sum over the momenta has been replaced by an integration over the first 
Brillouin zone of the square lattice. 

With the relation -£-r = coth(x/2) — 1, one obtains the following expression for 
the magnetization 

(S z ) = [- dkj dkyCotH^—)}- 1 . (130) 

TT Z JO JO 2 

This equation must be iterated to self-consistency in (S z ), which can then be com- 
pared with QMC (see Section 4.1.5). 

4.1.2. The Callen decoupling 

In this section, we discuss an attempt of Callen to improve the RPA. We do this 
because it is the basis of an approximate decoupling of the terms stemming from 
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the single-ion anisotropy (see Section 4.2.1). This generalisation of the Tyablikov 
(RPA) decoupling results from the ansatz 



((S?S+; Sr)) ~ (Sf)((S+; Sj)) - a (S; S+)((Sf; Sj)), (131) 

with a = t|^-; a — > corresponds to the Tyablikov (RPA) decoupling. Inserting 
this expression into the equation of motion and applying the spectral theorem leads 
again to a single-pole expression for the Green's function with a modified dispersion 
relation. The spectral theorem yields 



IX 2 Jo 



Callen 
k 



-1 



( SZ ) = L dk * L dk v coth(^ ) , (132) 

with 



ry rir rir T. O, .Callen 

^ allen = B+(S z )(j -J k )(l + - dkj d^^ooth(^— )). (133) 

71 JO JO Jq Z 



Again, eqn ()132j) must be iterated to self-consistency in (S z ). Although it takes 
some higher-order correlations are into account, the Callen approach is worse than 
RPA for the present case but still much better than a mean field (MFT) result (see 
Section 4.1.5). 

oo 

Derivation of the Callen dispersion relation (|133j) . 

In order to make the Callen decoupling plausible, consider two equivalent formulas for 
spin S = 1/2 

QZ Q Q— Q+ 

St = \(SfSr-SrSt). (134) 

Multiplying the first equation by a and the second by (1 — a), one can write the Green's 
function ((SfS^;S~})m the following form: 



((SfS+;Sr)) = aS((S+;Sr)) + 1(1 - a) ((S+S-S+; S~)) - i(l + a)((SrS+S+; S~)). 

(135) 

Now factorize the Green's functions on the right hand side as follows: 



((SfSrS+;ST)) ~ (SfSr)((S+; Sj)) + (SfS+)((Sr- SJ)) + (SrS+)((Sf; 5J», 
((SrSfS+;Sj)) * (S.St}((S+;Sr)} + (S.S+)((St;Sr}} + (S+S+}((Sr;Sr)). 

(136) 

Approximating the terms non-diagonal in the z-component of the spin by (S^~S^~) ~ 
and using the relation (S^S^) = 2(Sf) + (S~ S^), we obtain 

((S!S+-,Sr))^(Sf)((S+-,Sr))- a (SrS+)((St-,Sr)), (137) 
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which is expression (|131|) . Taking a = ^U- interpolates between the case a = 1, where 
the first of the equations (|134|) should be used for the decoupling at low temperatures 
({S z } ~ S), and a = 0, where the second formula should be used ({S z ) ~ 0). For arbitrary 
spins, arguments in favor of a = (which includes the spin S = 1/2 case) are given in 

El- 

Introducing the decoupling (|137|) into the equation of motion (|121|) yields 

( w _ B - (S*) Ju - a£ J« <SfS+)) Go -( w ) 
i i 

+ ((S Z ) Ju +aJ2Jil(SrS?))G lj (co) = 25^). (138) 
i i 

A Fourier transform to momentum space leads to 

u — B — (S Z )(J - J k ) - a± £(J q - J q+k )(5-5+) q ]G k (a;) = 2(S Z ), (139) 

where the Green's function is given by 

U - LOy. 

with 

fallen = B+ ( ^ )(Jq _ Jk) + a j_ £ (Jq _ J q+k) ^ 5 +> q . (141) 

The spectral theorem then determines 

{S ~ S+)k = 27 So J-oo ^^T (Gk(tJ + i5) - Gk{uJ - iS) = 2{SZ)q>k - (M2) 



with 



I -i / a. .C alien \ 

__ = i( coth( ^_)_l). (143) 



It remains to simplify the term proportional to a in the dispersion relation: 

= 2a(5 z )^ ]T J^l - e ik ( R ^^ )l £ e iq(Ri-R^ q 

= 2a<S*>(J -J k )^£^* q , ( 144 ) 
where we have made use of 

iVV q N^zJ 2 ^ q iVWo q ' 1 ' 

q q <»j> q 

and z is the number of nearest neighbours. This completes the proof for the Callen 
dispersion relation. 

oo 
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In Ref. jT7], Callen also derives a closed form expression for the magnetization 
for general spin S from the solution of a differential equation. The result is 

, on _ (S - $ k )(l + $ k )( 2S+1 ) + (S + 1 + $ k )4 25+1) n ^ 

{ } ( l + $ k) (25 + l)_^ + l) > ^ 



a formula which was also found by Pravecki [18J. For the treatment of general spin 
S, see also Refs. [T9] . 

4.1.3. Mean field theory (MFT) 

In mean field theory (MFT), which is frequently used as the simplest approxima- 
tion, one neglects correlations which lead to collective excitations (magnons). The 
essential approximation consists in writing operator products as 

S?Sf = (5? - - (£/» + S?(Sf) + (S?)Sf - (Sf><^> 

*S?(S?) + (S?)Sf-(S?)(S?) t (147) 



where the mean field assumptions S" — (iSf) and Sj ~ (Sf) have been made. 
Neglecting transverse expectation values ((Sf) = 0) as well leads to the mean field 
Hamiltonian 

h mft = Ju{sl)st - B J2S?+l £<^> (Sf), (148) 

kl I Z kl 

where the last term, being a constant, does not influence the equations of motion 
but has to be taken into account when calculating the intrinsic energy. 

In Green's function theory, the Hamiltonian H MFT leads without further ap- 
proximations to the equations of motion 

{to — B — (S z ) £ Jf^Gij = 2{S z )5 ij , (149) 

k 

whose Fourier transform to momentum space is 

(u-B-(S z )J )G k {u) = 2{S z ), (150) 

where Jq = zJ (z is the number of nearest neighbours; z = 4 for a square lattice) 
and 

G,(u) = , (151) 

with a momentum-independent dispersion relation 

co MFT = B + (S z )J . (152) 
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Because there is no momentum dependence in this relation, the /c-integration in the 
spectral theorem is trivial, so that 

- - - (S-S+) -1-iy 2{SZ) - - - 2{SZ) (153) 

from which 

i a MFT 

<^> = -tanh(^— ). (154) 

This result is obtained from the RPA result by setting Jk to zero in eqn ()126j) . 
thereby neglecting the k-dependence of the lattice. In MFT, it is only the number 
of nearest neighbours z that count. This is also true for the more complicated cases 
discussed later. The neglect of the k-dependence (this corresponds to the neglect of 
magnons) makes MFT much worse than RPA, as seen in Fig. ^of Section 4.1.5 . 

Because MFT is easily applied, often with qualitatively reasonable results, we 
quote a few papers where MFT is extensively used: in Refs. [20] and and refer- 
ences therein, the spin reorientation transition is treated and effective (temperature- 
dependent) lattice anisotropy coeefhcients are calculated; in Ref. [7%] . coupled ferro- 
antiferromagnetic layers are treated. 

4.1.4. The Mermin- Wagner theorem 



Mermin and Wagner [22] have shown quite generally that the pure Heisenberg model 
(without magnetic field and anisotropies) in less than 3 dimensions does not exhibit 
collective order at finite temperatures (for (S z ) — > the Curie temperature goes to 
0: T Curic -> 0). 

From the expressions derived above, we can see that RPA obeys the theorem 
whereas MFT violates it. Expanding the RPA expression for the magnetization 
(I13(jp for small (S z ) and B = 0, one obtains an expression for the Curie temperature 

r RPA - 1 ^ 1 > n n ^) 

1 Curie - v 2 ~ 2 p 11 rvr r> 2 ^ ptJL fir 1 ^ U ' y 100 ) 
j _j k 2 7T 2 JO ^xJo «Sj -J k JO "K-zJo k 2 +k 2 

This is so because the integral diverges at the lower boundary, which can be seen by 
expanding the square lattice dependence of Jk, eqn ()127j) . for small momenta. This 
means that the Mermin- Wagner theorem is obeyed in RPA. In three dimensions, 
the RPA expression gives a finite value for the Curie temperature and is used in 
ab initio calculations of the Heisenberg exchange interaction to determine the Curie 
temperature, see e.g. [23] . 

Calculating the Curie temperature from the MFT result (|154|) for B = by 
expanding for small (S z ) gives a finite Curie temperature 

ZcST = \j = \zJ = J, (156) 
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where we have taken z — 4, the number of nearest neighbours for a square lattice. 
This is in clear violation of the theorem. 



4.1.5. Comparing with Quantum Monte Carlo calculations 

In this section, we compare the temperature dependence of the magnetization of a 
ferromagnetic spin 5 = 1/2 Heisenberg monolayer on a square lattice obtained with 
the approximations just described with quantum Monte Carlo (QMC) calculations 
[To] , which are 'exact' within their statistical errors. The results are shown in Fig^ 
The RPA of Section 4.1.1 is the best and is a fairly good approximation. Before 
QMC calculations were available, it was not possible to check the quality of RPA. 
Although there are additional correlations taken into account in the Callen approach 
of Section 4.1.2, its results for spin S = 1/2 are not as good as those from simple 
RPA. In Ref. [17J, Callen argues that his decoupling should give better results 
for larger spin values, but there are no QMC calculations available to support his 
statements. The mean field theory (MFT) of Section 4.1.3 yields by far the worst 
results. This results from not taking collective excitations (magnons) into account, 
which is also the reason for the violation of the Mermin- Wagner theorem. 

We mention that RPA gives still better results for the magnetization when higher- 
order Green's functions with vertex corrections for their decoupling are included 
For quantities with transverse correlations, like the intrinsic energy or the specific 
heat, one has to go beyond RPA. See Section 5, in particular Ref. [SB] . 

4.1.6. The effective (temperature dependent) single-ion lattice anisotropy 

Lattice anisotropy coefficients are defined in an expansion of the free energy in 
powers of cos# [21], where 9 is the polar angle between the magnetization (S) and 
the normal to the film plane 



The anisotropy coefficients can be calculated by thermodynamic perturbation the- 
ory, where the Hamiltonian H = H + V is separated into an unperturbed part 
H consisting of the exchange coupling and the magnetic field and a perturba- 
tion V n = —K n J2i(Sf) n (n=2,4). Within first order perturbation theory, effective 
anisotropy coefficients can be defined as 



where the temperature dependence is introduced by the functions f n {T) which are 
expressed in terms of expectation values ((S z ) n )o for the unperturbed Hamiltonian 



F(T, 6) =F Q {T) - K 2 (T) cos 6 - K 4 (T) cos 4 6 - B ■ (S) . 



(157) 



K n (T) = K n f n (T) , 



(158) 




(159) 
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square monolayer S=1 /2 




Figure 1: The temperature dependence of the magnetization of a ferromagnetic 
Heisenberg monolayer for a square lattice with spin S = 1/2. Comparison between 
the 'exact' quantum Monte Carlo (QMC) result ^3] and the results obtained with 
MFT, RPA and Callen decoupling. We have used (a) J/B=20, (b) J/B=10 and (c) 
J/B=4 (from reference |14j). 

U(T) = |[35((S 2 ) 4 ) - (30S(S +l)-25)((^) 2 ) + 3S(S + 1)(S(S + 1)- 2)] . 

The moments are calculated with RPA and MFT in Ref. ^3] and the resulting 
temperature dependent coefficients are shown in Fig. (j2J) for S = 2 and S = 10, 
where K n = 1 and a scaling J — ► J/S(S + 1) and B — > S / 5 has been used. The 
resulting behaviour of the K n (T) calculated by RPA differs markedly from that 
obtained by MFT particularly at low temperatures: whereas the JC n (T) obtained 
with MFT show an exponential decay in this temperature range, those calculated 
from RPA decrease more rapidly and exhibit a nearly linear behaviour. The JC(T) 
calculated with RPA exhibit a much weaker dependence on the spin S than those 
calculated with MFT. 

In Section 4.2.2 we show that it is better to calculate the effective anisotropy 
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coefficients non-perturbatively by minimizing the free energy with repect to the re- 
orientation angle. 





T/e 



Figure 2: The temperature dependence of the effective lattice anisotropy coefficients 
K,2(T) and K/^T) of a square Heisenberg monolayer calculated with thermodynamic 
perturbation theory for MFT and RPA. We have used J/B = 100 and (a) S — 2 
and (b) S = 10. To allow for comparison between different spin values, we used the 
scaling J -> J/S(S + 1) and B -> B/S. 
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4.2. Ferromagnetic Heisenberg films with anisotropies for 
general spin S 

An isotropic Heisenberg model in less than three dimensions does not show spon- 
taneous magnetization at finite temperature, as explained by the Mermin- Wagner 
theorem [22] • Such an idealized system does not, however, exist in nature, since even 
the smallest anisotropy leads to a finite magnetization. This can be caused by an 
external magnetic field (as shown in the previous chapter), single-ion anisotropies, 
exchange anisotropies, or the magnetic dipole-dipole interaction. 

Many applications of GF-theory deal only with the magnetization in one direc- 
tion of space. They treat multi-layers but not all use the full power of the eigen- 
vector method outlined in Section 3.3. We mention only a few. Diep-The-Huang et 
al. |2H1 treat ferro- and antiferromagnetic multilayers but, instead of the eigenvec- 
tor method, they use Kramers rule for calculating the GF's. Schiller and Nolting 
j2Hl treat sc(100) and fcc(lOO) ferromagnetic Heisenberg spins with S = 7/2 using 
RPA for the exchange interaction and the Lines decoupling j2H] for the single-ion 
anisotropy. C. Cucci et al. j2Zj consider fee (100), (110) and (111) ferromagnetic 
Heisenberg films using RPA, the Lines decoupling and the eigenvector method. 

In the following, we do not restrict the magnetization to be in one direction 
of space because we are interested in the reorientation of the magnetization as a 
function of the temperature and film thickness. Therefore, we deal from the outset 
with a multi-dimensional case; the orientation of the magnetization in one direction 
and the monolayer then occur naturally as special cases. An essential complication 
connected with the reorientation problem is the occurrence of zero eigenvalues of the 
equation-of-motion matrix, which can be handled with the techniques developed in 
Section 3. 

We do not discuss papers dealing with the magnetic reorientation on the basis 
of a boson expansion, as e.g. Refs. [3U1 EI] who start with a Holstein-Primakoff 
transformation in lowest order, because the validitity of a linearized spin wave theory 
is limited to low temperatures only. 

4.2.1. The Hamiltonian and the decoupling procedures 

We consider a spin Hamiltonian consisting of an isotropic Heisenberg exchange in- 
teraction with strength Jjy between nearest neighbour lattice sites, an exchange 
anisotropy with strength D^, a second-order single-ion lattice anisotropy with strength 
i^2,fc, a magnetic dipole coupling with strength g^i and an external magnetic field 
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B = (B x ,By,B z 

n= - 



Z <kl> Z <kl> k 



I) 2 



- Y.{\B-St + \B + S- k +B z Sl 

k 

+ iT.^l^Sf + SiSf) -3(S*r„)(S,r„)). 

A kl r kl 



(160) 



Here the notation = S k ± iS k and B ± = B x ± %B V is introduced, where k and 
I are lattice site indices and < kl > indicates summation over nearest neighbours 
only. 

In order to treat the spin reorientation transition for general spin S, we need the 
following Green's functions: 



(161) 



where a = (+, — , z) takes care of all directions in space, 77 = ±1 refers to the 
anti-commutator or commutator Green's functions respectively, and n > 1, m > 
(m + n < 2S + 1) are positive integers. We follow the formalism of Section 3 by 
evaluating all formulas for the Hamiltonian ([160)1 . 
The exact equations of motion 



"G%™(u) = A%™ + (([S?, W]_; (S*) m (S~r)U 
with the inhomogeneities 



A 



{[s?,(s*) m (s-T] v ) 



(162) 



(163) 



are given explicitly by 



ujG 



±,mn 



±,mn 



A 



TE4(«; (S]r(Sj) n )) - ((S z k Sf; (S^ m (S. 

k 

±J2D ik ((S* k St;(S*) m (S-T)) 

k 

±K 2 Msfsi + s*sfy, (s;r(s;) n )} 



ULr ij,ri 



a z,mn 



J2M((SrS+ - S k St); (S*) m (S.) n )) 



2 W7 



+ 2^ ^,17 



(164) 



For the moment, we leave out the terms due to the dipole-dipole interaction, which 
we include later. 
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Once these equations are solved, the components of the magnetization can be 
determined from the Green's functions via the spectral theorem. A closed system of 
equations results from decoupling the higher-order Green's functions on the right- 
hand sides. For the exchange interaction and exchange anisotropy terms, we use a 
generalized Tyablikov- (or RPA-) decoupling: 

{{S a S fi. {S *T{SJY)) V ~ (S?)G%™ + (SfrGflZ". (165) 

The terms stemming from the single-ion anisotropy have to be decoupled differently, 
because RPA decoupling leads to unphysical results; e.g. for spin 5 = 1/2, the terms 
due to the single-ion anisotropy do not vanish in RPA as they should do because, in 
this case, ^2,i{(Sf) 2 ) is a constant and does not influence the equations of motion. 
In the appendix of Ref . [2E] , we investigate different decoupling schemes proposed 
in the literature, e.g. those of Lines or that of Anderson and Callen [32] • These 
should be reasonable for single-ion anisotropics small compared to the exchange 
interaction. We found the Anderson- Callen decoupling to be most adequate in our 
context. It treats the diagonal terms as they occur from the single-ion anisotropy in 
the same way that Callen J7| used in his attempt to improve the RPA. Consider eqn 
()131J) for i — I: add the term for ((S-Sf; ...)) and do the same for the corresponding 
expressions for G~' mn . Using S^Sf = S(S + 1) =F S? — S?Sf, one obtains 

(((s?sf + s;s?y 1 (sj) m (s;) n )) v 

This term vanishes for S = 1/2 as it should. 

In Section 4.2.5, we shall demonstrate a procedure for treating the single-ion 
anisotropy exactly by going to higher-order Green's functions. With this, single-ion 
anisotropies with arbitrary strength can be treated. This procedure is, however, 
tedious to apply for spins S > 1, whereas there is no problem when staying at the 
level of the lowest-order Green's function as discussed in the present section. 

Applying the decouplings ()165|) and ()166|) and a Fourier transform to momentum 
space, one obtains, for a ferromagnetic film with N layers, 3N equations of motion 
which can be written in compact matrix notation as 

(ul - T)G m " = A mn . (167) 

Qmn j g a sjV-dijxiensiona] Green's function vector and 1 is the 3N x 3N unit ma- 
trix. The Green's functions and the inhomogeneity vectors each consist of N three- 
dimensional subvectors which are characterized by the indices i and j, which are 
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now layer indices: 



f Gp mn (k,cu) \ I A^ mn \ 

A 

i AZ,mn j 



, 4 ( k >^) / 



A ran 



-.van 
I z,mn 



(168) 



The equations of motion are then expressed in terms of these layer vectors and 
tatrice: 

^ T\\ fn . . . Ti at 



3x3 submatrices T^- of the 3N x 3iV matrix T 



ul - 



r 12 

"21 



\ r^i t N2 ■ ■ ■ r NN j 













G 2j 




A 2j 5 2j 











1,...,N 



(169) 

The r matrix reduces to a band matrix with zeros in the sub-matrices, when 
j > i + 1 and j < i — 1. The diagonal sub- matrices Tu are of size 3x3 and have 
the form 

( Hf -H+ \ 

o -Hf h; 





V 



2 » 



(170) 



where 
W - 

Hf - 
H± = 



B* + (S!)(ju{q - 7k) + Dug) + (J M+1 + A,w)(^i) + ( J M-i + A,i-l)(^ti) 
1 



+^2,2(^)1 



2S 2 



[5(5+1) -(5^ 



5f + (S*) J«(g - 7 k) + Ji,i+i{St+i) + Ji,i-i(Sti) 
H? - (St) Dm . 



(171) 



For a square lattice and a lattice constant taken to be unity, 7k = 2(cos k x + cos k y ) 
and q = 4 is the number of nearest neighbours. 

If the dipole-dipole coupling is small compared to the exchange interaction, it can 
be treated in the mean field approximation (see e.g. the appendix of jH] and Section 
7.4 of this review). In this case, the dipole coupling leads to a renormalization of 
the magnetic field and one finds 



N 

Bf = B^ + Y,9iASf)T^\ 

N 

B * = fl«-2 5><,<S;)Tl i -'l; 

7 = 1 



(172) 



i.e. there is an enhancement of the transverse fields and a reduction of the field 
perpendicular to the film plane. 
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The lattice sums for a two-dimensional square lattice are given by (n = \i 

I 2 - n 2 

rpn _ \ ~> « " 

~ (I 2 + m 2 + n 2 ) 5 / 2 ' 



-J\) 
(173) 



The indices Im run over all sites of the j layer excluding terms with I +m +n = 0. 
The 3x3 non-diagonal sub-matrices for j — i ± 1 are of the form 



J^Sf) 

\ ) —\Jij{Sf) 







(174) 



Now the system of equations of motion is completely specified. 

The case (iQj 7^ 0, = 0) has been treated in Ref. j2H] for a monolayer and 
in Ref. [5| for the multilayer by the eigenvector method. In this case H i = H i . 
In the case ((ify = 0, D^- 7^ 0), treated in |37|, one has Elf 7^ ii^ , which leads to 
additional dependencies on the momentum vector k, requiring a refinement of the 
treatment. We discuss these cases separately in the following subsections. 



4.2.2. Approximate treatment of the single-ion anisotropy 



For the single-ion anisotropy, one can use eqn ()49|) directly because the term R°L° 
turns out to be independent of the momentum k. The +, —,z components of the 
vector C k n are, however, not independent, i.e. there are not enough equations to 
solve for the unknowns. The remedy is to supplement eqn (}4*9"|) with the regularity 
conditions (j43j) . 

For illustration, consider the monolayer. For = and K 2 7^ and H ± = H ± , 
the eigenvalues of the equation-of-motion matrix T ()170|) and eigenvector matrices 
R and L can be determined analytically. The eigenvalues are uq = 0, d± = ±i?k = 
±\ / H + H^ + H Z H Z . The right eigenvectors are arranged so that the columns 1,2 
and 3 correspond to the eigenvalues 0, +E k and — E k respectively: 



R 



/ 


H+ 


-(E k +H*) 




\ 




H z 


H~ 


H~ 






H~ 


(Ew-H*) 


-{E^+H*) 






H z 


H+ 


H+ 




V 


1 


1 


1 


/ 



(175) 



and the left eigenvectors are arranged in rows 1,2,3 corresponding to the eigenvalues 
0, +E k , —E k : 



( 



AEl 



2H~H Z 
(E k + H Z )H~ 
\ (E k — H Z )H~ 



2H+H Z AH Z H Z \ 
(E k - H Z )H + 2H~H+ 
-(E k + H Z )H+ 2H-H+ j 



(176) 
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With the knowledge of L° the regularity conditions (J45j) are 

/ A+J ran \ 

L°A_i = = (H~H Z , H + H z , 2H Z H Z ) AZ'" m ■ 
((S-)J(q- lk ) + B-)(S z ) 



(177) 



For m = 0, n = 1, with Atf 1 = 2(S Z ), AZf = and Azf 



H 

IT" 



(S z ) 



B z + (S z )J(q - 7k ) + K 2 2(S*)(l - ^[S(S + 1) - (S*S')] 



Solving for (S ) and taking the complex conjugate, 

B ± 



(S z ). 



(178) 



(179) 



B* + K 2 2(S*){1 - ^[S(S + 1) - (S'S*)^ 

Thus, once (S z ) and (S Z S Z ) have been calculated, the transverse correlations follow 
from the rgularity condition. Note that the prefactor jp- does not depend on the 
momentum vector k. 

The lowest spin for which the single-ion anisotropy has an effect is S = 1 (for 
S =1/2 the anisotropy term is a constant and does not contribute to the equations 
of motion). In this case, only equations of motion with (m = 0, n = 1) and (m = 
l,n = 1) are needed to determine the correlations (S z ) and (S Z S Z ). The regularity 
conditions with m + ?2<2S'+l = 3 suffice to express all remaining correlations as 
functions of (S z ) and (S Z S Z ) (for more details see Appendix B of Ref. 28j). 

For the monolayer with general spin S, comparison of ()178j) with ()179|) shows 
that 

(180) 



El 

Z ' 



with Z = B z + K 2 2(S Z )[1 - 5^(5(5 + 1) - (S Z S Z ))). The regularity conditions 
()177|) can therefore be written for general m, n in the form 



2ZA z,mn = A +,rnn B - + J^ran^ 



:isr 



Using the ^-component of equation ()47j) for the monolayer, one obtains a relation 
between the correlations in momentum space 



2 I L^_l(s z ) m {S-) n S z ) 
Z 2 

_ 1 i+,mn ^k B~ rff k 

~ 2 ~ l H z Z \-H z 



B' 



coth( 



■((S z ) m (S-) n S + ) 



+ -AZ? 
2 - 1 



~ ~Y 

E k B+ r E 



((S z ) m (S-) n S-) (182) 



H z Z 



^+coth(- 2 



Note that all correlation functions in this equation are written in a standard form 
where powers of S z are written to the left of the powers of S~: 



C{m,n) = ((S z ) m (S-) n ). 



(183) 



40 



The relations [S z , (S~) n ]- = -n(S-) n and S-S + = S(S + 1) - S z - (S z ) 2 allow us 
to express all correlations in terms of the C(m,n): 

((S z ) m (S-) n S z ) = nC(m,n) + C(m+l,n) , 

({S z ) m {S-) n S + ) = (S(S + l)-n{n-lj)C(m,n-l)-(2n-l)C{m + l,n-l) 

-C(m + 2,n-l) , 

((S z ) m (S~) n S-) = C{m,n + 1). (184) 

The inhomogeneities can also be expressed in terms of the C (m, n) using binomial 
series: 

A_li = —nC (m, n) , 

A +,mn = _ q m _ ( S z ) m ) S~ S + + 2S*{S Z - l) m + (n - l)(n + 2S z ){S z ) m ] (S^)"" 1 ) 

m / m \ 

= S{S + l)Y,\ j (-iyC(m-i,n-l) + (2n + m)C(m + l,n-l) 
m+1 ( m + 1 \ 

+ E . (-l)' +1 C(m + 2-z,r2-l) + n(n-l)C(m,n-l) , 



i=2 



AZT n = ([(S z + l) m -(S z ) m }(S~) n+l ) = J2 y . JC(m-i,n + l) . (185) 

The regularity conditions for all m and n can be written in terms of correlations 
defined in the standard form by inserting equation ()185j) into equation (j!81|) : 

r m ( m \ 

2ZnC{m,n) = B-[S{S + l)Y,[ . I (-l)*C(m - i, n - 1) 



m+l 



m + 1 



+ (2n + m)C(m + l,n-l)+ ^ I . ] (-l) i+1 C(m + 2 - i,n - 1) 



m 



+n(n-l)C{m,n-l) +B + Y,\ ]C(m-i,n + l) . (186) 



For a given spin 5 and given values of C(m, 0) for m < 2S + 1 this set of linear 
equations is solved for C(m,n > 0) for all m + n < 2S + 1. The resulting values 
are then checked for consistency by insertion into the 25* equations ()182|) using 
equations ()184|) and ()185jl . The solution consists of moments ((S z ) p ) (p—l,...,2S) 
for which eqns ()182j) are self-consistently fulfilled. Note that the highest moment 
{(S z ) 2S+l ) has been eliminated in favour of the lower ones through the relation 
U Ms (S z -M s ) = 0. 

Multilayers are treated simply by adorning the correlations and the quantity Z 
with a layer index i: 

Z% = B? + Ji,i+i(S? +1 ) + J i) i_i(5'f_ 1 ) 
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+K 2>i 2(S*)(l - ^[S(S + 1) - (S*S?)}) 



;i87) 



An alternative method of solution of the present problem is first to eliminate 
the null-space by a singular value decomposition (SVD) of the equation-of- motion 
matrix T (|17U|) and then using equation (J77)l directly. The advantage is to reduce 
the dimension of the problem by the number of zero eigenvalues. 

The monolayer for S — 1 is suitable [SB] for demonstrating the procedure because 
the SVD of r and the vector v of eqn (|77jl can be obtained analytically. In order 
to consider a reorientation of the magnetization in the xz-plane, we put B y = 0, so 
that H = H x . The T-matrix can be expressed as a product of three matrices: 

r = UWV = uwv. (188) 



Proceeding as in Section 3.5, it is a simple exercise to obtain the three factors: 
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Now everything is specified and one can solve equation (|77jl. There is the technical 
problem that the vectors at neighbouring fc-values in general have arbitrary phases 
but this can be overcome with the help of the smoothing procedure described in 
Section 3.5. 

In Fig. El we show as an example the results for the spin S = 1 Heisenberg 
monolayer. It does not matter whether one uses eqn (J4*9~|) or eqn (|77jl : the results 
are the same, as they should be. 

This is not the case in later examples (see Sections 4.2.3 and 4.3.2) where it 
is necessary to use the singular value decomposition to deal with a momentum- 
dependent factor R°L°. 
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Heisenberg monolayer: J=1 00, K 2 =1 , B x =0.4 
field-induced reorientation of the magnetization 




For spin S > 1 and multilayer systems, one can use both methods as long as 
R°L° is independent of the momentum, but the singular value decomposition and 
the matrix v now have to be determined numerically 

In the next figures, we show further examples from Ref. PJ. Fig. 0] shows 
normalized magnetizations (S z )/S and (S x )/S for a monolayer as functions of the 
temperature for all integral and half-integral spin values between S = 1 and S = 6 
calculated with Green's function theory. The reorientation temperature Tjf depends 
slightly on S. The inset shows the corresponding results in mean field theory for 
spins S = 1,2,7/2, and 11/2. In this case, the reorientation temperature does not 
depend on S. Note the very different temperature scale, which is due to the missing 
magnon correlations in MFT. 
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Figure 4: Normalized magnetizations (S z ) / S and (S x )/S for a monolayer as func- 
tions of the temperature for all integral and half-integral spin values between 5=1 
and S = 6 calculated with Green's function theory. The reorientation temperature 
Tf depends slightly on S. The inset shows the corresponding results in mean field 
theory for spins S = 1,2,7/2, and 11/2. The reorientation temperature for MFT 
does not depend on S. 




Fig. shows the equilibrium reorientation angle as a function of the temperature 
for the systems of Fig. 0] calculated with GFT. The inset shows the corresponding 
results for MFT. 
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Figure 5: Equilibrium reorientation angle as a function of the temperature for the 
systems of Fig. HI calculated with GFT. The inset shows the corresponding results 
for MFT. 
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Fig. El shows the sublayer magnetizations (S?) as functions of the temperature 
for thin ferromagnetic films with N layers and spin S = 1. The reorientation 
temperature Tr for the different films can be read off from the curve in the N — T 
plane, where (Sf) = 0. 




Figure 6: Sublayer magnetizations (S?) as functions of the temperature for thin 
ferromagnetic films with N layers and spin S — 1. The reorientation temperature 
Tr for the different films can be read off from the curve in the N — T plane, where 

(Si) = o. 
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In Fig. the average equilibrium reorientation angle Go is shown as a function 
of the temperature for different film thicknesses. N is the number of layers in each 
film and T% are the reorientation temperatures at 0q = 90°. 
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Figure 7: The average equilibrium reorientation angle Go as a function of the tem- 
perature for different film thicknesses. N is the number of layers in each film and 
T$ are the reorientation temperatures at 6 = 90°. 



If one is interested in the effective (temperature-dependent) lattice anisotropy 
coefficient ^(T), one should not use the thermodynamic perturbation theory dis- 
cussed in Section 4.1.6, but rather a non-perturbative approach in which the free 
energy is minimized with respect to the layer-dependent reorientation angles df 
dFi(T)/d9i\ e . = 0, where 9 i0 are the equilibrium reorientation angles. The effective 
anisotropy of a film consisting of iV layers is 



N 



(192) 
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cos6 0ji (B x + J M+1 M m (T) sin0 o ,i+i + J^iM^ mn0 o ,i-i + Tf n ) 
- sin 6 ,i(B z + Ji, t+ iM l+l (T) cos0 O) i+i + J M -iM M _! cos0 Oii -i - 27™ s 



Here M»(T) = J (Sf) 2 + (Sf) 2 and 0o,i — arctan((5'f ) / (/S'f )) are determined from 
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Figure 8: Average effective anisotropy TCziT) / '£2(0) as a function of the temperature 
and film thickness N. The inset demonstrates the different functional dependence on 
T for layers with N=l and N=19 if the temperature is scaled with the reorientation 
temperature. 

the magnetization components and 

N 

rr yin.^i,^^,/ 1 '' 1 , (194) 

where T'* - - 7 ' are dipole lattice sums, see (|173|) . 

In Fig. |Sl average effective anisotropies ^(iV, T)/K,2(N, 0) of films with different 
thicknesses N are shown as functions of the temperature. 

Up to now, we have used S + , S~ , S z as the basic operators to define the Green's 
functions suitable for treating the reorientation of the magnetization in the case 
of uniaxial anisotropies. When there are anisotropies in all directions of space, 
it is more natural to start with the operators S x ,S y ,S z , because this treats the 
three directions of space on an equal footing. This is done in Ref. where 
the Anderson-Callen decoupling of the single-ion anisotropy terms is invoked for all 
directions of space. A formal advantage is that the equation-of-motion matrix turns 
out to be hermitean. Generalising a formula due to Callen ^7] leads to analytical 
expressions for the first and second moments of the spin operators. Reorientation 
transitions and effective (temperature-dependent) anisotropies are calculated for 
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various 3D and 2D cases. 

The GF theory is used in Ref. [31] to investigate the interplay between a single- 
ion easy-plane anisotropy and the dipole-dipole interaction for a Heisenberg mono- 
layer with the Hamiltonian (K 2 > 0) 

h = -jJ2 + K 2 E(^) 2 - B E s i + Hdipole 

<ij> i i 

The Tyablikov decoupling is used for the exchange and the dipole-dipole inter- 
actions and the Anderson- Callen decoupling for the single-ion anisotropy. An inter- 
esting result is that the easy-plane anisotropy alone cannot stabilize the long-range 
ferromagnetic order at finite temperatures, one needs the dipole-dipole interaction 
in addition in order to do so. 

The spin reorientation problem is also investigated in Refs. [331 In these 
papers, a single-ion anisotropy is used and the dipole-dipole interaction is approx- 
imated by the dipole demagnetization energy. The exchange interaction and de- 
magnetization energy terms are treated by the Tyablikov (RPA) decoupling and a 
decoupling due to Lines [2H| is applied to the single-ion anisotropy terms. However, 
instead of calculating the longitudinal and transverse components of the magnetiza- 
tion vector, only the z-component of the magnetization is calculated as a function 
of the temperature. For the multi-layer case, the vanishing of the gap in the corre- 
sponding spin-wave spectrum at a particular temperature is interpreted as the onset 
of the reorientation transition. Effective (temperature-dependent) anisotropies are 
also calculated within this approximation. 
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4.2.3. Treating the exchange anisotropy 



Although the formalism described in Section 4.2.1 looks very similar for the single- 
ion and exchange anisotropy, the direct application of the standard spectral theorem 
is not possible because the term R°L° in eqn (|47)l turns out to be momentum 
dependent owing to the fact that for the exchange anisotropy H i ^ Hf in eqn (j!71|) . 
Thus the Fourier transform in the second term of eqn (|4"7j) cannot be performed. 

In Ref. {37] we found by intuition a transformation which eliminates one of the 
rows of R°L° in the equation 



C = R 1 £ 1 L 1 A_ 1 + R°L°C 



(195) 



thus allowing the corresponding row to serve as an integral equation of the eigen- 
vector method. 

The transformation is found to be 
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with T L T = 1. 

Applying this transformation to equation ()195|) 

T _1 C = T~ 1 R 1 f 1 L 1 TT _1 A_ 1 + T^RoLqTT^C 

and inserting the analytical eigenvectors R and L for the monolayer 
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(e k -H z )H x 2H X H X , (199) 
-(e k + H Z )H X 2H X H X j 

the second component of the vector T" 1 R L TT _1 C transforms to zero and one 
obtains, together with the regularity conditions (|45|). the integral equations for the 
correlations for each (m, n) pair. 

The eigenvector method immediately generalizes to the case of N layers if the 
transformation T is extended to 3iV dimensions by constructing 3iV x 3iV-matrices 
with sub-matrices ()196|) on the diagonal. 

The intuited transformation (|196|) can also be found systematically enlisting the 
help of the singular value decomposition of the T-matrix as described in Sections 
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3.5 and 4.2.2. This automatically yields some momentum-independent components 
of a row vector v 3 - which enables the Fourier transformation (|78|). 

This procedure also works for the case of coupled ferro- and antiferromagnetic 
layers described in Section 4.3.2. 

Some typical results for systems with exchange anisotropy are shown in Figs, 
and I n Fig- EH the magnetization (S z ) and its second moment (S Z S Z ) are 
plotted as functions of the temperature for a ferromagnetic spin 5=1 Heisenberg 
monolayer for a square lattice with an exchange interaction strength of J = 100 and 
an exchange anisotropy strength of D = 0.7. It is interesting to note that there is 
practically no difference in the magnetization curves when using an Anderson- Callen 
decoupled single-ion anisotropy, once its strength is fitted to an appropriate value, 
K 2 = 1.0. This makes it difficult to decide which kind of anisotropy is acting in an 
actual film. 
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Figure 9: The magnetization (S z ) and its second moment (S Z S Z ) of a ferromagnetic 
spin S — 1 Heisenberg monolayer for a square lattice as functions of the temperature, 
comparing a GFT calculation using an exchange anisotropy (D = 0, open circles) 
with a single-ion anisotropy (K 2 = 1.0, solid dots). The corresponding MFT results 
are also shown. Note the different Curie temperatures. 
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A novel feature occurs with the introduction of the magnetic dipole coupling: the 
eigenvalues and eigenvectors of the r-matrix become complex above a certain tem- 
perature, i.e. below a certain value of (S z ). This behaviour occurs quite naturally 
in the theory. It has nothing to do with a damping mechanism and has to be taken 
seriously in order to obtain the results of Fig. Because the r-matrix is real, its 
eigenvalues and eigenvectors, if complex, occur pairwise as complex conjugates, and 
the integral equations to be solved must be real. 

In Fig. the components of the magnetization (S z ) and (S x ) and the absolute 
value S for a fixed magnetic field B x = 0.3 are shown as functions of the temperature 
for a ferromagnetic spin S = 1 Heisenberg monolayer for a square lattice. Also 
shown are the equilibrium reorientation angle, 0o and the critical reorientation 
temperature, Tr, at which the in-plane orientation is reached. The small horizontal 
arrow indicates the value of (S z ) below which complex eigenvalues occur. The results 



reorientation as function of the temperature 
at fixed B*=0.3 J=1 00 g=0.066 D=0.7 S=1 




T R temperature 

Figure 10: The components of the magnetization (S z ) and (S x ) and its absolute 
value S for a fixed magnetic field B x = 0.3 as functions of the temperature for a 
ferromagnetic spin S = 1 Heisenberg monolayer for a square lattice. The exchange 
interaction strength is J = 100, the exchange anisotropy strength is D = 0.7 and the 
strength of the magnetic dipole coupling is g = 0.066, a value corresponding to Co. 
Also shown are the equilibrium reorientation angle, Go, and the critical reorientation 
temperature, Tr, at which the in-plane orientation is reached. The small horizontal 
arrow indicates the value of (S z ) below which complex eigenvalues occur. 

for the exchange anisotropy and the single-ion anisotropy for spins S > 1 and multi- 
layers look very similar, as seen by comparing references jU] with jSTj. We therefore 
do not show the corresponding figures here. 
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Similar results for general spin S using the exchange interaction together with the 
dipole coupling in the mean field approximation are found in Ref. jHHj, but for the 
monolayer only. We also mention Ref. jUJ], where the competition of the exchange 
anisotropy with the single-ion anisotropy is investigated for a ferromagnetic S=l 
Heisenberg monolayer, and Ref. jH], where the formalism is applied to multilayers 
and higher spin values with parameters that are different at the surface and the 
interior of the film. In Ref. [12] the spin reorientation transition for a ferromagnetic 
Heisenberg monolayer with exchange interaction, exchange anisotropy and dipole 
dipole interaction is treated with the RPA, where, however, a somewhat artificial 
temperature dependence of the exchange interaction had to be used in order to 
obtain a favorable comparison with experiment. 
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4.2.4. Susceptibilities 

In reference [75] Jensen et al. report measurements of the parallel and transverse 
susceptibilities of a bi-layer Cobalt film having an in-plane uniaxial anisotropy. They 
analyse their results with the help of a many-body Green's function theory assuming 
an exchange anisotropy and a value for the spin of S — 1/2. Here, we generalize 
their theoretical model, extending it to multilayers and arbitrary spin. We discuss 
not only the exchange exchange anisotropy but also the single-ion anisotropy. A 
comparison of the two cases allows an evaluation of the robustness of the theoretical 
conclusions as well as possibly identifying any qualitative differences which might 
enable an experiment to discern which type of anisotropy is acting in a real film. Ac- 
cordingly, we investigate the parallel and transverse susceptibilities for arbitrary spin 
in multilayer systems. In keeping with the earlier work [75J , we use the Green's func- 
tion formalism and neglect the dipole-dipole interaction, since it is nearly isotropic 
for the in-plane case. The theory is formulated in complete analogy to Sections 4.2.2 
and 4.2.3 (where an out-of-plane magnetization was discussed), the only difference 
being that the applied magnetic fields allow only an in-plane magnetization. For a 
detailed description, refer to references [TH] and [77]. 

The adequate decouplings for the in-plane situation are the same as for the 
out-of-plane case: whereas a RPA decoupling is reasonable for the terms coming 
from the exchange interaction and the exchange anisotropy, it leads to incorrect 
expressions for the single-ion anisotropy terms. For the latter we therefore use the 
method proposed by Anderson and Callen [32] at the level of lowest order in the 
Green's function hierarchy. This is certainly an adequate approximation for small 
anisotropies, as we have shown in Ref. [IB] for the case of an out-of-plane single- 
ion anisotropy of a monolayer by comparing with 'exact' Quantum Monte Carlo 
calculations. In addition to Sections 4.2.2 and 4.2.3, we refer the reader to the 
literature for a discussion of the roles of the single- ion- [U1I2E1E] and exchange |37j- 
anisotropies with respect to reorientation of the magnetization of a ferromagnetic 
film with an out-of-plane anisotropy as a function of temperature and film thickness. 

The susceptibilities with respect to the easy (Xzz) and hard (xxx) axes are cal- 
culated as differential quotients 

X ZZ = ((S z (B z ))-(S z (0)))/B z 

Xxx = ((S*(B*))-(S*(0)))/B*, (200) 

where a value B z ^ = 0.01/ S turns out to be small enough. We compare numer- 
ical results obtained with the single-ion anisotropy with those from the exchange 
anisotropy. As the single- ion anisotropy is not appropriate for S = 1/2, we show 
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Figure 11: The magnetization (S z )/S of a ferromagnetic spin S = 1 Heisenberg 
monolayer for a square lattice shown as a function of the temperature. Green's 
function (indicated by RPA) calculations with exchange anisotropy D = 5 (crosses) 
and with single-ion anisotropy (K 2 = 5.625), (open squares) in the Anderson- Callen 
approximation are compared. Also shown are the quantities 100 * (S X )/(S + 1) for 
the exchange anisotropy and 100 * (S x )/S for the single-ion anisotropy; the factor 
100 is introduced to make the curves visible. The corresponding results for mean 
field (MFT) calculations are also displayed. 




results for S > 1. In an attempt to obtain universal curves (i.e. independent of the 
spin quantum number S), we scale the parameters (B x ( z \ J, D) in the Hamiltonian 
as B x ^/S = B x{ - Z \ J/S(S + 1) = J and D/S(S + 1) = D (D being the strength 
of the exchange anisotropy). We also scale the strength of the single- ion anisotropy 
according to K 2 /{S - 1/2) = K 2 . 

In order to compare results obtained with the single-ion anisotropy with those 
of the exchange anisotropy, we set the strength of the single-ion anisotropy to 
K 2 = 5.625 for a square lattice monolayer with spin S=l, so that the easy axis 
magnetization (S z )/S lies as close as possible to the magnetization obtained with 
the exchange anisotropy (D = 5) used in [76 . The exchange interaction parameter 
is J = 100 and there is a small magnetic field in the x-direction, B x = 0.01/ S, which 
stabilizes the calculation. The comparison is shown in Fig. ^2 It is surprising that 
the results for the easy axis magnetization (S z ) are very similar over the whole tem- 
perature range although the physical origin for the anisotropies is very different. An 
analogous result was observed for the out-of plane situation discussed in Ref. |37| . 
For the exchange anisotropy, the hard axis magnetization is a constant below the 
Curie temperature, whereas for the single-ion anisotropy, it rises slightly up to the 
Curie temperature. In Ref. [TH], it is shown analytically that the hard axis magne- 
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RPA+MFT: <S I >/S, 1 00.<S">/(S+1 ) 
monolayer: S=1 12 ,1 ,3/2,2,3,4,6,1 3/2 




Figure 12: The magnetizations {S*)/S of spin 5 = 1/2, 1,3/2,2,3,4,6, 13/2 Heisen- 
berg monolayers for a square lattice as functions of the temperature, from Ref. [To] . 
Results from Green's function (RPA) calculations are compared with results from 
mean field theory (MFT) using the exchange anisotropy strength D = 5. Also shown 
is the hard axis magnetization, which scales to a universal curve 100 * (S X )/(S + 1), 
where the factor 100 is introduced to make the curves visible. 

tization for the exchange anisotropy is universal for a scaling (S X )/(S + 1). For the 
single-ion anisotropy, a scaling (S x )/S is found to be more appropriate. Comparison 
with the corresponding mean field (MFT) calculations, obtained by neglecting the 
momentum dependence of the lattice, shows the well-known shift to larger Curie 
temperatures (by a factor of about two for the monolayer with the present choice of 
the parameters) owing to the omission of magnon excitations. 

In Figs. and EH we show the easy and hard axes magnetizations for a mono- 
layer with different spins S. Whereas in Fig. H*2lone observes a nearly perfect scal- 
ing for the RPA calculations with the exchange anisotropy (S = 1/2,1,3/2,2,3, 
4,6,13/2) and a universal Curie temperature Tc{S) for RPA and MFT, this is 
not the case for the corresponding results with the single-ion anisotropy shown for 
S = 1, 3/2, 4, 5 in Fig. although the violation of scaling is not dramatic. 

Turning to the inverse easy and hard axes susceptibilities Xzz an d Xxli we find 
very similar results for the exchange anisotropy and the single-ion anisotropy. In 
particular, in the paramagnetic region (T > T Curie ), the inverse susceptibilities as 
a function of temperature are curved owing to the presence of spin waves, whereas 
the corresponding MFT calculations show a Curie- Weiss (linear in the temperature) 
behaviour. There is slightly less universal behaviour for the single-ion anisotropy 
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Figure 13: The magnetization {S z )/S of ferromagnetic spin S = 1,2,3/2,5 Heisen- 
berg monolayers for a square lattice as a function of the temperature for Green's 
function (RPA) calculations using the single-ion anisotropy strength of K 2 = 5.625 
and the corresponding results of mean field theory (MFT). Also shown are the quan- 
tities 100 * (S x )/S; the factor 100 is introduced to make the curves visible. 



(Figs. ITU and rH^) than for the exchange anisotropy (Figs. 2 and 3 of Ref. [23 )• This 
is connected with the fact that the exchange anisotropy exhibits universal values for 
the Curie temperatures T^ F (S) and T^ FT {S), which is not strictly the case for 
the single-ion anisotropy, (Fig. fl3*)l . We were also able to show analytically in Ref. 
[76] that Xxx * S(S + 1) 1S universal for T < Tc for the exchange anisotropy; this is 
not the case for the single-ion anisotropy. The only difference is in the curves for the 
imperfectly scaled Green's function results for x7z- f° r ^ ne exchange anisotropy, the 
curve with the lowest spin value lies to the left of the curves with the higher spin 
values, whereas the converse is true for the exchange anisotropy. This is not a very 
pronounced effect and does not lead to a significant difference between the results 
for the various anisotropies. In treating multilayers with the exchange anisotropy in 
Ref. [7H], we considered only the case S = 1/2. The single- ion anisotropy term in 
the Hamiltonian is a constant for 5 = 1/2; therefore it is not active when calculating 
the magnetization, so we have to use a larger spin here. In the following, we use spin 
S=l as an example but we also have results for S > 1 which scale with respect to 
the spin in the same way as in the monolayer case. The Curie temperatures for the 
multilayers N = 2, 19 (for N=19 one is already close to the bulk limit) are only 
slightly lower for the single-ion anisotropy than those calculated for the exchange 
anisotropy. 

Some results are shown in Figs. EH an d El 111 order to avoid cluttering the 
figures, we restrict ourselves to a multilayer with N=9 layers and spin S — 1. For 
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Figure 14: 'Universal' inverse easy axis susceptibilities Xzz * S(S + 1) °f an i n ~ 
plane anisotropic ferromagnetic square lattice Heisenberg monolayer as functions of 
the temperature for single-ion anisotropy and spins S = 5,2,3/2, 1. Compared are 
Green's function (RPA) and mean field (MFT) calculations. 



N > 9 the corresponding curves would shift only slightly in accordance with the 
saturation of T c with increasing film thickness. We display only the RPA results for 
the multilayer (N=9) and compare with the RPA monolayer (N=l) result. Again, 
there is no significant difference in the results for both anisotropies. We do not 
plot the corresponding mean field results, which are shifted to higher temperatures 
and, in the paramagnetic region, show only (a linear in T) Curie- Weiss behaviour, 
whereas the RPA results have curved shapes owing to the influence of spin waves, 
which are completely absent in MFT. 

Although both kinds of anisotropies are of very different physical origin, it is pos- 
sible, by fitting the strengths of the anisotropies properly, to obtain nearly identical 
values for the easy axis magnetizations over the complete temperature range for a 
spin 5 = 1 monolayer. Using the parameters obtained in this way for monolayers 
with higher spin values and for multilayers, we looked for differences in the results 
of calculations with both kinds of anisotropies. 

By using scaled variables we find a fairly universal behaviour (independent of 
the spin quantum number S) of easy and hard axes magnetizations and inverse 
susceptibilities. Universality holds better for the exchange anisotropy; e.g. we find 
a universal Curie temperature T C (S) for RPA and MFT. The scaling is not as 
perfect for the single-ion anisotropy, but there are no dramatic deviations which 
might enable an experiment to distinguish between the two types of anisotropies. It 
is sufficient to do a calculation for a particular S and then to apply scaling to obtain 
the results for other spin values. In principle the measurement of the hard axis 
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Figure 15: 'Universal' inverse hard axis susceptibilities Xxx * + 1) °f an i n ~ 
plane anisotropic ferromagnetic square lattice Heisenberg monolayer as functions of 
the temperature for single-ion anisotropy and spins S = 5,2,3/2, 1. Compared are 
Green's function (RPA) and mean field (MFT) calculations. 



susceptibility together with the Curie temperature allows one to obtain information 
about the parameters of the model, the exchange interaction and the anisotropy 
strengths. One should, however, keep in mind that the quantitative results of the 
present calculations correspond to a square lattice. They could change significantly 
for other lattice types. Further changes could result from the use of layer- dependent 
exchange interactions and anisotropies. Such calculations are possible, because the 
numerical program is written in such a way that layer-dependent coupling constants 
can be used. 

A general result is that there are no qualitative differences for the calculated 
observables (easy and hard axes magnetizations and susceptibilities) between the 
single-ion anisotropy on the one hand and the exchange anisotropy on the other 
hand. Therefore, it is not possible for us to propose an experiment that could 
decide which kind of anisotropy is acting in a real ferromagnetic film. 

We mention also a paper by Yablonskyi |43j . who derives analytical expressions 
for the static susceptibility and for correlation functions for ferromagnetic and an- 
tiferromagnetic Heisenberg monolayers with general spin (no anisotropies) on the 
basis of the Tyablikov (RPA) decoupling. 
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Figure 16: The inverse easy axis susceptibilities Xzz °f ferromagnetic films in RPA 
for spin S = 1 for a monolayer (N=l) and a multilayer (N=9) as functions of the 
temperature for single-ion and exchange anisotropies. 
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Figure 17: The inverse hard axis susceptibilities Xxx °f ferromagnetic films in RPA 
for spin 5=1 for a monolayer (N=l) and a multilayer (N=9) as functions of the 
temperature for single-ion and exchange anisotropies. 



4.2.5. Exact treatment of the single-ion anisotropy 

Up to now we have worked at the level of the lowest-order GF's, where approximate 
decoupling schemes lead to closed systems of integral equations which are solved 
self-consistently. In this subsection, we show that a closed system for the terms 
stemming from the single-ion anisotropy is attainable without any decoupling by 
going to higher-order GF's 03], generalizing the work of Devlin (35] • By taking 



60 



advantage of relations between products of spin operators |17] , one can show that 
the hierarchy of the equations of motion is automatically closed with respect to 
the anisotropy terms . In this way, an exact treatment of the single-ion anisotropy 
results, i.e. anisotropies of arbitrary strength can be treated, whereas e.g. the 
Anderson-Callen deccoupling of the second-order single-ion anisotropy terms is only 
reasonable for anisotropies small compared to the exchange interaction. The terms 
due to the exchange interaction must still be decoupled by a generalized RPA. 

We develop the general formulation for a spin-Hamiltonian consisting of an 
isotropic Heisenberg exchange interaction between nearest neighbour lattice sites, 
Jm, second-order and fourth-order single-ion lattice anisotropies with strengths K 2 , k 
and K^k respectively, a magnetic dipole coupling with strength g k i and an external 
magnetic field B = (B x ,By,B z ): 

ft = _ 9 E Jki(SkSf + s^s") 

Z <kl> 

-E^(^) 2 -E^4, fc (^) 4 
k k 

- E (W S k + \ B+S k + b z si) 

k 

+ lT, g 4{ r li(SkS^ + S z k Sn-S(S k -r kl )(S l -r kl )) , (201) 

Z kl r kl 

where Sf = Sf ± iSf and B ± = B x ± iBy, k and I being lattice site indices and 
(kl) indicates summation over nearest neighbours only. We have added a fourth- 
order anisotropy term for which we had no decoupling procedure available when 
working at the level of the lowest-order GF's. 

To allow as general a formulation as possible (with an eye to a future study of 
the reorientation of the magnetization), we formulate the equations of motion for 
the Green's functions for all spatial directions: 

= ((St;S?)) u 
G^(u) = ((Sr;S?)) u (202) 

Instead of decoupling the corresponding equations of motion at this stage, as we 
did in our previous work [2E1II3 5 we add equations for the next higher-order Green's 
functions: 

Gtp^(u) = ((S z St + StS?)) = (((2S?-l)St;Sj))^ 

G ij z,T ( u ) = ((s;sr+srsi)) = ((sr(2s;-i) ] s?)) u 
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((StSt;SJ}) w (203) 

(((6S!Si-2S(S+l));S?)) u . 

The particular form for the operators used in the definition of the Green's functions 
in Eqs. (I203j) is dictated by expressions coming from the anisotropy terms. Ter- 
minating the hierarchy of the equations of motion at this level results in an exact 
treatment of the anisotropy terms for spin S = 1 , since the hierarchy for these terms 
breaks off at this stage, as will be shown. The exchange interaction terms, however, 
still have to be decoupled, which we do with RPA-like decouplings. 

For the treatment of arbitrary spin S, it is necessary to use 4S(S + 1) Green's 
functions in order to obtain an automatic break-off of the equations-of-motion hier- 
archy coming from the anisotropy terms. These are functions of the type Gfj with 
a = (z) n (+) m and a = (—) m (z) n , where, for a particular spin S, all combinations of 
m and n satisfying (n + m) = 2S must be taken into account. There are no Green's 
functions having mixed + and — indices because these can be eliminated by the 
relation S^S ± = S{S + 1)^5*- {S z ) 2 . 

Here we treat only the spin S = 1 monolayer, for which there are 8 exact equa- 
tions of motion for the Green's functions defined in (J202j) and (|203|) . 

The crucial point now is that the anisotropy terms in these equations can be 
simplified by using formulae which reduce products of spin operators by one order. 
Such relations are derived in Ref. |T7j : 

2S-m 

(S-) m (S z ) 2S+1 ' m = (S-) m 6$ s ' m) (S z y, 

i=0 

2S-m 

(s z ) 2S+1 ~ m (s + ) m = sl s ' m) {s z y{s + ) m . (204) 

i=0 

The coefficients 8\ S are tabulated in Ref. f|7J for general spin. For spin 5 = 1, 
only the coefficients with m = 0,1,2 occur: 8^'°^ = 8^'°^ = 0; 8^ = 1,8q = 
0,8^ = 1,8^ = 1. 

These relations effect the reduction of the relevant Green's functions coming from 
the anisotropy terms in the equations of motion 
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The higher Green's functions coming from the anisotropy terms are thus expressed 
in terms of lower-order functions already present in the hierarchy; i.e. with respect 
to the anisotropy terms, a closed system of equations of motion results, so that no 
decoupling of these terms is necessary. In other words, the anisotropy is treated 
exactly. For higher spins, S > 1, one can proceed analogously. 

No such procedure is available for the exchange interaction terms, which still 
have to be decoupled. For spin S — 1, an RPA-like approximations effects the 
decoupling: 



(206) 



Note that we have constructed the decoupling so as not to break correlations having 
equal indices, since the corresponding operators form the algebra characterizing the 
group for a spin 5 = 1 system. For spin 5 = 1, this decoupling model leads to 8 
diagonal correlations for each layer i: 

(Si~), (S i ), (Sf), (S^S^), (S i S i ), (S£Si~), (Si Sf), (SI SI). 

These are determined by the 8 decoupled equations. Performing in addition a two- 
dimensional Fourier transformation to momentum space results in a set of equations 
of motion which, in compact matrix notation, are 



(ul — r)G T = A T , 



(207) 



where and are 8-dimensional vectors with components G a,T and A a ' T and 

a = +, — , z, z+, — z, ++, ,zz; 1 is the unit matrix. The 8x8 non-symmetric 

matrix T is 
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H+ 


-2H Z 
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((25 z - 1)5+) J k 
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(s-s-)j k 
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-H z 







-((2S Z - l)S+)J k 
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-H+ 





2H Z 
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-2(5-5->J fc 





H- 
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3(5- (25- -l))J k 


-3((2S* - l)S+)J k 





-3H- 


3H+ 






with the abbreviations 



Jk 

K 2 



= B« + (S a )J(q- lk ), 



B a + (S a )Jq, 
K 2 + K 4 . 



a = +, 



a = +, -, z 



(209) 
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For a square lattice with a lattice constant taken to be unity, 7k = 2 (cos k x + cos k y ), 
and q = 4 is the number of nearest neighbours. For spin S = 1 and S = 3/2, 
the K± term in the Hamiltonian leads only to a renormalization of the second- 
order anisotropy coefficient: K 2 (S = 1) = K 2 + K 4 and K 2 (S = 3/2) = K 2 + \K± 
respectively. Only in the case of higher spins, S > 2, are there non-trivial corrections 
due to the fourth- order anisotropy coefficient. 

If the theory is formulated only in terms of G~, there is no equation for determin- 
ing the (S + S + ) occuring in the T— matrix. It is for this reason that we introduced 
G + in Eq. (J2(J2j) . for which the T— matrix is the same, so that, in general, one can 
take a linear combination of G + and G~ and their corresponding inhomogeneities: 



G = (1 
A = (1 

Hence, the equations of motion are 

(ul - 



a)G~ +aG~ 
a)A~ + aA" 1 



(210) 



r)G = A, (211) 
= (1 — a)C~ + aC + can be determined. The 



from which the desired correlations C 
parameter a is arbitrary (0 < a < 1). 

An examination of the characteristic equation of the T-matrix reveals that 2 of 
the eigenvalues are exactly zero, so that the term R°L° is needed when applying 
the eigenvector method of Section 3.3. The eigenvector method then yields for the 
correlations in configuration space (i = 1, ...,8): 

C * = h f dk * f*i(EE RljSj^L^At + RlL%d) . (212) 

17 J0 J0 1=1 j=lk=l j=l 

Without loss of generality, the field component B y can be set to zero, which leads 
to the symmetry requirements: (S + ) = (S~), (S + S + ) = (S~S~) and (S Z S + ) = 
(S~S Z ); i.e. there are only 5 independent variables defining 8 correlations C, i.e. 
the system of equations is overdetermined. This problem can be overcome with a 
singular value decomposition: define a vector consisting of the five relevant quantities 

( (S~) \ 
(S z ) 

(SS-) ■ (213) 

(s-s z ) 
V (s z s z ) ) 

Then, the correlations C can be expressed as 



C 



u: 



u c v 



(214) 
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The 8x5 matrix u c may be written in terms of its singular value decomposition: 

u c = UWV, (216) 

where W is the 5x5 diagonal matrix of singular values which here are all > for 
< a < 1. U is an 8x5 orthogonal matrix and V is a 5 x 5 orthogonal matrix. 
From Eqs. (j2T2jl and QZTty it follows that 



u c v = R 1 £ 1 L 1 A + R°L°(u c v + u°) - u°. 



(217) 



To get v from this equation, we need only multiply through by u c 1 = VW 1 U, 
which yields the system of coupled integral equations 



v = u, 1 (R 1 S 1 L 1 A + R°L°(u c v + u°) 



(218) 



or more explicitly with i = 1 , . . . , 5 

8 



Vi 



k=l 71 j=l 1=1 

2 5 8 

i=i P =i k=i 

This set of equations is not overdetermined (5 equations for 5 unknowns in the 
present example ) and is solved by the curve-following method described in Appendix 
B. 

As an example we investigate the magnetization as a function of the second- 
order anisotropy strength and the temperature for a spin S — 1 square monolayer, 
putting the dipole coupling and the magnetic field equal to zero. In this case the 
magnetization is in the ^-direction only, (S z ). The results are shown in Fig. [THl 
together with those from the Anderson-Callen decoupling. There is rather good 
agreement for small anisotropies, which, however, worsens as K 2 increases. An- 
other difference concerns the second moments, (S Z S Z ), which approach the value 
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K£ exact (open circles); decoupling {lines) 



T 




Figure 18: The spin S = 1 monolayer with exchange interaction strength J = 100. 
Comparison of GFT calculations for (S z ) and {S Z S Z ) as functions of the temperature 
for various anisotropies using the exact treatment of the anisotropy (open circles) 
and the Anderson- Callen decoupling of Section 4.2.1 (small dots). 

(S Z S Z )(T — ► T C urie) = 2/3 for the Anderson-Callen decoupling (see Ref. |28j), 
whereas in the exact treatment, the values of {S Z S Z )(T — > Tcurie) are larger than 
2/3. Estimates for the Curie temperature, as e.g. in Refs. jSHI or |1B] . give reason- 
able values only for small single-ion anisotropies. 

To show the difference between the new model and the Anderson-Callen decou- 
pling more clearly, we compare in Fig. EH the Curie temperatures obtained from 
MFT, the Green's function theory with the exact treatment of the anisotropy and 
the Green's function theory with the Anderson-Callen decoupling of Refs. [2E1 HI] • 
For small anisotropy, there is only a slight difference between the two GFT results 
which, in contrast to MFT, obey the Mermin- Wagner theorem. However, for large 
anisotropy, the GFT results deviate from one another significantly: for K2 — > 00, 
the Anderson-Callen result diverges, whereas the exact treatment approaches the 
MFT limit. This is shown analytically in the appendix of Ref. jH]. 

Unfortunately, we have not been able to solve the full reorientation problem with 
the exact treatment of the single-ion anisotropy with the tools developed in Section 
3.5, because of numerical difficulties. 

When using the Anderson Callen decoupling we obtained rather good results 
when the external field is in the direction of the anisotropy as long as the anisotropy 
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Curie temperatures: MFT and GFT 
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Figure 19: Comparison of the Curie temperatures calculated with the exact treat- 
ment of the anisotropy, the Anderson-Callen decoupling and MFT. The first two 
approaches fulfill the Mermin- Wagner theorem: Tc — > for — > 0, whereas MFT 
does not. For large anisotropies, the exact treatment approaches slowly the MFT 
result (as also can be shown analytically [44J), whereas the Anderson-Callen decou- 
pling leads to a diverging T c 



is small enough (K 2 < 0.1J). This is seen by comparing with Quantum Monte Carlo 
calculations [IB] • The approximation is much worse when the field is applied perpen- 
dicular to the anisotropy. A considerable simplification and an improvement of the 
results concerning the reorientation is reported in Ref. [JO], where the Anderson- 
Callen decoupling is made in a frame which is rotated with respect to the original 
one and in which the magnetization is in the direction of the new z-axis. The reori- 
entation angle is determined from the condition that the magnetization commutes 
with the Hamiltonian in the rotated frame. In this connection see also Ref. |5D] . 
who also apply the approximate Anderson-Callen decoupling in a rotated frame. In 
Section 4.4.1 we treat the spin reorientation with an exact treatment of the single-ion 
anisotropy by working also in the rotated frame. 
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4.2.6. The importance of spin waves in the Co/Cu/Ni trilayer 



The importance of spin waves can be demonstrated in Co/Cu/Ni trilayers, where 
two magnetic layers are separated by a non- magnetic spacer layer. In an experiment, 
the magnetization of Ni in a Ni/Cu bilayer and a Ni/Cu/Co trilayer is measured as 
a function of the temperature [53] . Figure I2H shows a shift to higher temperatures 
of the magnetization curve of Ni for the trilayer system (dots) as compared to the 
Ni magnetization in the bilayer system (crosses). This shift is largest at the Curie 
temperature. In the figure, results from Green's function theory are also shown. 
A Heisenberg exchange interaction and a dipole-dipole interaction can explain the 
observed shift with realistic strengths [54] for the interlayer coupling 0.5 < Ji n t er < 
3.0 (53J, assuming an in-plane magnetization, whereas MFT (owing to the neglect of 
spin waves) needs unrealistic strong values for Jinter- For more recent experimental 
results concerning Co/Cu/Ni/Cu(100) layers and a comparison with GFT, see Refs. 
[53} 157)] and references therein. 
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Figure 20: The measured shift of the Ni magnetization curve for the trilayer 
Co/Cu/Ni system (dots) as compared to the Ni magnetization curve for the bi- 
layer Ni/Cu system (crosses). Green's function theory (lines) can explain this shift 
with realistic strengths for the interlayer coupling Ji n t er [SB]; whereas MFT needs 
unrealistic strong values for Ji nte r- 
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4.2.7. Temperature dependence of the interlayer coupling 

The interlayer coupling between the ferromagnetic layers of Section 4.2.7 is caused 
by the spin-dependent reflection of spacer electrons at the magnetic/non-magnetic 
interface leading to a spin-dependent interference and to a renormalisation of the 
density of states and the free energy within the non-magnetic spacer. The coupling 
may then be ferromagnetic or antiferromagnetic, oscillating with respect to the 
spacer thickness with a period depending on the Fermi surface of the spacer. The 
amplitude and the phase of the coupling is determined mainly by the spacer thickness 
but may also be influenced by the interface roughness, disorder etc.. 

The sources of the temperature dependence of the interlayer coupling are dis- 
cussed in Refs. |5T| I58|: 

(i) A part of the temperature dependence is induced by the smearing out of the 
Fermi surface of the spacer, as proposed in j^I]. ('spacer effect') 

(ii) The temperature dependence also stems from altering the properties of the 
magnetic layers through spin wave excitations [52] ('magnetic layer effect'), which 
can affect the interlayer coupling. 

In Ref. inn], ferromagnetic resonance (FMR) experiments lead to an effective 
Ji(T) ~ T 3 / 2 dependence. Both mechanisms contribute and the dominant mech- 
anism cannot be deduced directly. In Ref. |58j . an alternative analysis of FMR 
measurements is proposed that could distinguish between both mechanisms. 
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4.3. Antiferromagnetic and coupled ferromagnetic- antifer- 
romagnetic Heisenberg films. 

A Green's function theory of antiferromagnetic (AF) and coupled ferro- and anti- 
ferromagnetic (AF-AFM)-films relies on periodic structures and therefore requires 
the introduction of sublattices in which periodicity is guaranteed. We start with a 
description of an antiferromagnetic monolayer in subsection 4.3.1 and follow this in 
subsection 4.3.2 with a general formulation in terms of sublattices, which allows a 
unified treatment of FM , AFM and FM-AFM multilayer-systems. 

4.3.1. The antiferromagnetic spin S = 1/2 Heisenberg monolayer. 

According to the Mermin- Wagner theorem [22] the two-dimensional antiferromag- 
netic or ferromagnetic Heisenberg monolayers with exchange interaction alone can- 
not show a finite magnetization. In order to obtain a finite magnetization for the 
antiferromagnet, one can either introduce an artificial (staggered) field with op- 
posite directions for the up and down spin sublattices [3] or one can introduce 
anisotropies. We use an exchange anisotropy and demonstrate how the magneti- 
zation of an antiferromagnet can be calculated with many-body Green's function 
theory. The essential step is the introduction of separate sublattices for the up and 
down spins. 

Consider the Hamiltonian 

n = -\Y. MS k St + S z k S?) -W D* kl S z k S?, (220) 

Z <kl> Z <kl> 

where the exchange interaction and the exchange anisotropy strengths are negative 
(J kl < and D z kl < 0). 

We only consider the magnetization in z-direction. The equation of motion for 
the relevant Green's function in energy space 

G±- = ((St;Sr)) (221) 

is 

uG±- = 2<^)<% + (([S+,H]; Sj)). (222) 

Again, we adopt the Tyablikov (RPA)-decoupling of the higher-order Green's func- 
tions occurring on the right-hand side: 

((SfS+;Sr))^(Sf)((S+;Sr)) . (223) 
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This leads to the equation 

(u; - £(J a + D* U )(S?))G±- + (Sf) E JuGtf = Z(Sf)5 ir (224) 
i i 

We now introduce sublattice indices (m, n) for the up (u) and down (d) spins. 
Four equations of motion corresponding to the pairs (i„,j m ) = (u,u), (d,u), (u,d) 
and (d, d) result. 

Fourier transforms to momentum space for the sublattices each consisting of N/2 
lattice sites are 



JV k 

— Ve-*^-^) = 5- ■ 
JV k 

|^ e i(k-k R im=4k/; (225) 

where the subscripts of the Green's functions in momentum space G mn (k) now 
denote sublattice indices and not lattice sites. 

Because {S z )d = —(S z ) u for an antiferro magnet, the 4 equations of motion decou- 
ple to two identical pairs of equations which determine (S z ) u or (S z ) d respectively. 
Before replacing (S z ) d by —(S z ) u , the equations for G+~(k) and G d ~(k) are 

(a; - (S'U J uu (0) - J uu (k) + D* uu (0)) - (S*) d ( J ud (0) + £&(0))G+"(k) 

+ (5*> u J tl<l (k)G+-(k) = 2(5*> fl 

(u;-(S%(J du (0) + D z u (0))-(S%(J dd (0) - Mk) + D z dd (0))Gj-(k) 

+ (S z ) d J du (k)Gt~(k)= 0. (226) 

Restricting the coupling to nearest neighbours only implies that all interaction terms 
with equal sublattice indices are zero: J uu = D z uu = J dd = D z dd = 0. After replacing 
(S z )d by —(S z ) u , the matrix equation is 

oo + (S z ) u (J ud (0) + D z ud (0)) (S%J ud (k) \ ( G+"(k) \ = ( 2(S% \ 

-(S%J ud (k) co - (S%(J ud (0) + D* ud (0)) ) { G+T(k) ) { ) 

(227) 

For a square lattice with lattice constant a — 1, 

J ud (k) = J k = ^^J iuld e- ik ^- n ^ = 2 J(cos k x + cos k y ) 
J ud (0) + D z ud (0) = J Z = A(J + D Z ). (228) 
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Eliminating G^ u (k) from the two equations yields 

2(S z ) u (u-(S z ) u J z ) 



Q-(k) 



(u + (S z ) u J§)(u - (S%J§) + (S>)IJZ 



(229) 



with the poles 



u; lj2 = ±(S') u y/((JS) 2 - Jl) (230) 

From the spectral theorem, after integrating over the first Brillouin zone and using 
the relation {S~ S + ) u = 1/2 — (S z ) u for spin S = 1/2, the following equation for the 
sublattice magnetization (S z ) u for the up-spins results: 



i r 

12 + — / , 

-a 1 Jo 



coth(/3o;i/2) = . 



(231) 



10 U)\ 

This must be iterated to self-consistency in (S z ) u . Results for J = —100 and 
D z = —0.1, 1.0, —10.0 are shown in figure I2T1 



Antiferromagnetic Monolayer with Spin=1 12 
Exchange anisotropy D z : MFT and RPA 
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Figure 21: The sublattice magnetization of an antiferromagnetic Heisenberg mono- 
layer with exchange anisotropy as a function of the temperature for RPA and mean 
field (MFT) calculations for the parameters J = —100 and D z = —0.1, —1.0, —10.0. 



For the RPA result, the value of the sublattice magnetization at zero temperature 
is well below its saturation value of {S z ) u = 1/2, contrary to the situation for the 
ferromagnet. This is due to quantum fluctuations. The mean field limit , obtained 
by setting Jk = 0, does not show this suppression and also contradicts the Mermin- 
Wagner theorem by having a finite Neel temperature for D z — > 0. This theorem 
(T^eei — > for D z — > 0) is obeyed in the RPA calculation as can be seen by deriving 
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the Neel temperature from equation (J231)) by taking first the limit (S z ) — > and 
then D z -> 0. 

An analytical approximation to the Neel temperature results from a partial frac- 
tion decomposition of the expression obtained after expanding the hyperbolic cotan- 
gent in eqn (J231|) for small sublattice magnetization. Then the remaining integrals 
are expanded around k x = k y = or around k x = k y = tt repectively with the result 

T N « - -7rJ 2 . (232) 

The values T N (D* = -0.1) = 36.9, T N (D Z = -1.0) = 50.6 and T N (D Z = -10.0) = 
80.1 are only slightly higher (less than 10%) than the results of the exact calculations 
shown in figure |2"T1 

The extension to AFM multilayers can be found in Refs. jHZl EH] • 
The two-dimensional spin S — 1/2 Heisenberg antiferromagnet for a square 
lattice with nearest neighbour exchange interaction and dipole-dipole coupling (no 
anisotropy) is treated by Pich and Schwabl in Ref. [HI], where they use linear spin 
wave theory by applying the Holstein-Primakoff transformation [U2]- They obtain 
better results [63J for the Neel temperature (i.e. closer to experimental data) when 
applying GFT along the lines of Callen [TZj. In later papers, they use the same 
formalism to treat two-dimensional honeycomb antiferromagnets [HI] and to study 
the influence of the dipolar interaction in quasi-one- dimensional antiferromagnets 
on a hexagonal lattice 
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4.3.2. A unified formulation for FM, AFM and FM-AFM multilayers 

In this section we treat the coupled FM-AFM system in detail jHS] > introducing sub- 
lattices for both the AFM and FM parts. It will then be self-evident that each part 
by itself can be described as a special case by choosing the signs of the parameters 
appropriately. This shows that FM, AFM and coupled AFM-FM systems can be 
handled uniformly within the same formulation. 

There is previous work in which Green's function theory treats the coupling 
of ferromagnetic layers to antiferromagnetic layers: in reference [UH], a bilayer is 
investigated and reference jTD] treats an extension to multilayers. In both cases, 
only a collinear magnetization is considered. In reference |71j . a ferromagnetic film is 
coupled to an antiferromagnetic layer; however, the orientation of the magnetization 
of the antiferromagnet is frozen. Other work considers an antiferromagnetic coupling 
between ferromagnetic layers [721 113 Ell • 

In our discussion here, we allow a non-collinear magnetization, where the reorien- 
tation of the magnetizations of the ferro- and antiferromagnetic layers is determined 
by the interlayer coupling as in the MFT approach of [78] . We restrict ourselves 
to Heisenberg systems with spin S = 1/2 with an exchange anisotropy. This is 
not an essential restriction: references [73 [77j show for ferromagnetic layers that 
through an appropriate choice of anisotropy parameters the exchange- and single- 
ion anisotropies yield very similar results and that an appropriate scaling leads to 
universal magnetization curves for different spin quantum numbers. Below, we ex- 
amine in detail the magnetic arrangement of the simplest system: a perfectly ordered 
bilayer consisting of a FM monolayer that is coupled to an AFM monolayer. 

The starting point is an XXZ-Heisenberg Hamiltonian consisting of an isotropic 
Heisenberg exchange interaction with strength Jjj between nearest neighbour lattice 
sites, exchange (non-localized) anisotropies in the x- or z-directions having strengths 
Dfj and Dfj respectively and an external magnetic field B = (B x , 0,B Z ) confined to 
the film plane, which is the xz-plane: 
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<i?> <«j> 



- Y.( BX Si+ BZ S Z k)- (233) 



Again, S'f = Sf±iSf and < zj > indicates summation over nearest neighbours only, 
where i and j are lattice site indices. Because there is no field perpendicular to the 
film plane (B y = 0), the reorientation of the magnetization can only occur in the 
xz-plane. For the FM-AFM bilayer we choose the anisotropy of the ferromagnetic 
layer in the ^-direction, Dfj, and the anisotropy for the antiferromagnetic layer in 
the x-direction, D%. 

For 5 = 1/2, the required commutator Green's functions are 

G$-(u>) = ((S? ] Sr)) u , (234) 

where a = (+, — , z) takes care of all directions in space. A generalization to spin 
quantum numbers S > 1/2 is effected in a straight-forward way by introducing 
G% mn = ((Sf; (S?) m (Sy) n )) with m + n < 2S + 1 (m > 0; n > 1; m,n integer) as 
in Section 4.2.1. 

The equations of motion for the Green's functions in the energy representation 

are 

wGg» = A«- + «[Sf , H]; SJ})^ (235) 
with the inhomogeneities 

, (236) 



where (...) = Tr(...e~ /3W )/Tr(e~ /3W ) denotes the thermodynamic expectation value. 

In order to obtain a closed system of equations, the higher-order Green's func- 
tions on the right hand sides are decoupled as in Section 4.2.1 by a generalized 
Tyablikov- (RPA) decoupling 

((S?S%; S-)) v ~ (S?)Gfc + (S%)G?r. (237) 

After introducing two sublattices per layer, the resulting equations are Fourier trans- 
formed to momentum space according to eqns. ()225|) . yielding 
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±(B z + J2(S;)(J mp (0) + D z mp (0))) G r. 



1±- 

P 

P Z 

z v 



T(B* + J2(S;)(J mp (0) + D x mp {0)))G^ n 

p 

p 

-\{B* + £<^>(J mp (0) + D x mp {Q)))G + mn 

p 

+ n(Sm) Jmp(k)Gp n 
Z P 

+\{b* + J2(s;)(j mp (o) + D* mp (o)))G-- 

p 

~(S^)T, J mp^)G-. (238) 
z p 

For a square lattice with lattice constant do = 1, one has four nearest-neighbour 
intralayer couplings with sublattice indices n, m from the same layer 

Jmniy) = QO Jmn i JmnQ*-) = 7o(k) Jmn i 

(239) 

D%(0) = q D% , £^(k) = 7o (k) D£* , 

with the intralayer coordination number g = 4 and the momentum-dependent 
Fourier factor 

7o(k) = 2 (cos k x + cos k z ) . (240) 

Correspondingly, for the nearest neighbour interlayer couplings, m and n now being 
sublattice indices from different layers, one has 

^mn(O) — Qmt -J'mt i <^mn(k) = 7int(k) J ln t , 

(241) 

For sc stacking, the interlayer coordination number and the corresponding Fourier 
factor are given by 

9int =7int(k) = 1 . (242) 

For fee or bec stacking, 

q mt = 4 and 7int (k) = 4 cos(k x /2) cos(k z /2). (243) 

The mean field approximation is obtained by neglecting the Fourier factors, i.e. 
7o(k) =7int(k) = 0. 
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By choosing the appropriate signs of the exchange interaction and the exchange 
anisotropy coupling constants, one can treat ferromagnetic, antiferromagnetic and 
mixed systems with coupled FM and AFM layers. 

The general formalism is valid for any number of layers and sublattices. If Z is 
the total number of sublattices of the system, the dimension of the set of equations 
()238j) is 3Z 2 . We restrict ourselves here to the investigation of the bilayer, so that 
there are four sublattices and the system of equations (|238jl is of dimension 48 with 
a corresponding Green's function vector. Closer inspection reveals that the system 
of equations has the following substructure 

/ / T \\ 




V 



o r 

o o 
V o o 



r o 
or; 



G 2 
G 3 

Vg 4 ; 



A 2 
A 3 
V A 4 J 



(244) 



where the diagonal blocks T are identical 12 x 12 matrices, whose explicit form can 
be read off from equations (|238j) . The sublattice Green's functions G n (n = 1, 2, 3, 4) 
are vectors of dimension 12 consisting of 4 subvectors, each of dimension 3: 

( G ln \ 

G2n 



G, 

G3„ 
\ G4 n J 

where the 3-component vectors are 

( G+- \ 



n 



1,2,3,4 



(245) 



G r 



C — 

mn 

G z ~ 



m 



1,2,3,4 . 



(246) 



The inhomogeneity vectors have the same structure: 
/ A ln 5 ln \ 



Ar 



' 2(5*) \ 


V -(s x m ) ) 



m, n 



1,2,3,4 . 



(247) 



A2 n 02n 
As n 03 n 
\ A 4n 54„ J 

The big equation (|244j) of dimension 48 for the bilayer can therefore be replaced by 
4 smaller equations of dimension 12: 



(ul -T)G n = A n for n =1,2, 3, 4 . 



(248) 



It turns out that the 12 x 12 T-matrix has 4 zero eigenvalues. In this case we can 
use the formalism of Section 3.5, where the singular value decomposition of T leads 
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to a system of integral equations for the correlations C n (k) corresponding to the 
GF's G n (see eqn (|77|)): 

= J dk(r£ 1 lvA_ 1) „ - vC n (k)) n = 1, 2, 3, 4. (249) 

Section 3.5 explains how to find a k-independent vector v having a layer structure, 
i.e. v = (0, .., 0, v n , 0, ., , 0). In this way, the non-diagonal correlations disappear 
from those rows in equation (J249|) corresponding to v n and the k-integration can be 
performed: / dkvC n (k)= v / dkC n (k) = vC n . In the present case v n is given by 

Vn = ((-^ > --L | l)tf lni (-L,--^,l)tf aB , 

(-±=, ~, l)5 3n , (-J=, --^=, l)5 4n ) , n = 1, 2, 3, 4. (250) 

Putting equation ()250|) into equation (j249|) yields 4 equations which contain the 
8 magnetization components implicitly. The necessary additional 4 equations are 
obtained from the regularity conditions (f79*j) 

J dkL A n = /" dkfi A n = , n = 1, 2, 3, 4, (251) 

which are obtained from the regular behaviour of the commutator Green's functions 
at the origin. The Uo are the eigenvectors of the singular value decomposition 
spanning the null-space of the matrix T. The resulting 8 integral equations are 
solved self-consistently by the curve-following method described in detail in the 
Appendix B. Note that the Uo are determined numerically only up to an orthogonal 
transformation. To ensure proper behaviour as a function of k, uo must be calibrated 
at each k. A procedure for doing this is indicated in Section 3.5 and presented in 
detail in an appendix of reference [T3j . 

We now present results for the bilayer ferromagnet, the bilayer antiferromagnet 
and the coupled ferro- and antiferromagnetic bilayer. All calculations are for an in- 
plane orientation of the spins of both layers. In each case we compare the results of 
Green's function theory (GFT) with those of mean field theory (MFT) obtained by 
putting the momentum-dependent terms equal to zero. In order to see the effects of 
the interlayer coupling most clearly, we use different exchange interaction strengths 
for each layer: 

(a) FM-FM: J 1FM = 100, J 2FM = 50, 

(b) AFM-AFM: J 1AFM = -100, J 2A fm = -50, 

(c) FM-AFM: J FM = 100, J AFM = -50. 

Because of the Mermin- Wagner theorem (221 > anisotropies are required in the Green's 
function description: we take D z = +1.0 for FM layers and D x = —1.0 for AFM 
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layers. These values are appropriate for 3d transition metal systems. For a compen- 
sated interface, the magnetizations of the FM and AFM layers are almost orthogonal 
to each other even at T = because of the interface exchange interaction J int . We 
choose the FM magnetization to be oriented in the ^-direction and the AFM mag- 
netization in the x-direction. Our particular choice of the anisotropies supports this 
arrangement not only at T = but also at finite temperatures. For other choices 
of anisotropies, the magnetic arrangement could be different. The interlayer cou- 
pling is assumed to be positive for the ferromagnetic bilayer and negative for the 
antiferromagnetic bilayer. For the coupled FM-AFM system, both signs are used. 
We consider three interlayer coupling constants with strength J- mt = 30, 75, 160, 
one smaller than the weakest exchange interaction, one larger than the strongest 
exchange interaction and one in between. 
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The ferromagnetic and the antiferromagnetic bilayers 

Results for the FM and AFM bilayers are presented in this subsection in order to 
have a basis for discussing the differences from the coupled FM-AFM bilayer. 




( c ) MFT: bi-layers FM-FM and AFM-AFM 

JlFM = "JlAFM = 1 00 J2FM = ~J2AFM = 50 




temperature 



Figure 22: (a) Green's function theory (GFT) for the ferromagnetic bilayer: The 
sublattice magnetizations are displayed as a function of the temperature for different 
interlayer couplings J int = 30 (dotted), 75 (dashed), 160 (solid). The exchange inter- 
action and anisotropy constants are Jifm = 100, J2FM = 50, D* FU = 1.0, -D| FM = 
1.0. 

(b) GFT for the antiferromagnetic bilayer: The sublattice magnetizations are 
displayed as a function of the temperature for two interlayer couplings J int = 
—30 (dotted), —160 (solid). The exchange interaction and anisotropy constants 

are Jiafm = —100, J2AFM = —50, -Diafm = ~~ 1-0, -D2AFM = — 1-0. 

(c) Mean field theory (MFT) for the ferromagnetic and antiferromagnetic bilay- 
ers with identical parameters: Ji( 2 )fm = |^i(2)Afm|, ^i(2)fm = |A(2)Afm|> 

<^intFM = "^intAFM ■ 



In figure I2*2*h. we show the sublattice magnetizations of the ferromagnetic bi- 
layer as a function of the temperature for three interlayer couplings calculated with 
Green's function theory (GFT). The magnetization profiles are different for the two 
layers (the magnetization is larger for the layer with the larger exchange interaction) 
but end in a common Curie temperature, which increases with the strength of the 
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interlayer coupling: T Curie = 50.66, 55.24, 60.04. 

For the antiferromagnetic bilayer, the parameters are the same as for the fer- 
romagnetic bilayer except for a sign change. In figure 122b we show the sublattice 
magnetizations of the antiferromagnetic bilayer for two interlayer coupling strengths 
calculated with Green's function theory. To avoid clutter, we have left out the result 
for the intermediate interlayer coupling strength. The corresponding magnetization 
curves lie in between those of the other couplings. At low temperatures one observes 
clearly the well-known reduction of the magnetization due to quantum fluctuations, 
which are missing in MFT, see figure 122b. Since |Jiafm| > |^2Afm| this reduction 
is larger for the first layer. With increasing temperature the magnetization curves 
of the two layers cross each other (a fact which was first observed by Diep [6*7] ) 
and finally end in a common Neel temperature. A larger interlayer coupling leads 
to a larger suppression of the magnetization at low temperatures and to a larger 
Neel temperature. Whereas with the present choice of parameters the magnetiza- 
tion profiles of the FM and AFM bilayers are rather different at low temperatures, 
the critical temperatures turn out to be identical: T Curie = T Nfel (cf. figures l22T a) 
and (b)), a fact already discussed by Lines [8131 . 

For comparison, we show in figure [22(c) the results of mean field theory (MFT) 
with the same parameters. The magnetization profiles as well as the critical tem- 
peratures are identical for the ferromagnetic and antiferromagnetic bilayers. As is 
well known, the Curie (Neel) temperatures (Tc ur ie(Neei) = 102.10, 107.25, 123.16 ) 
are much larger (with the present choice of parameters by about a factor of 2) in 
MFT owing to the missing magnon excitations. In MFT the Curie temperature 
is not very sensitive to the anisotropics as long as they are much smaller than the 
exchange interaction. In GFT, however, the sensitivity is very much greater because 
of the Mermin- Wagner theorem [22] (Tc U rie(Neei) — > for D z ^ x > — > 0). Also, the effect 
of the interlayer coupling on the magnetization profiles is much stronger in MFT 
than in GFT. 
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The coupled ferro-antiferromagnetic bilayer 

This is the most interesting case. We consider here two in-plane magnetization corn- 





Figure 23: (a) Green's function theory(GFT): the sublattice magnetizations of the 
ferro- and antiferromagnetic sublattices are displayed as a function of the tempera- 
ture for different interlayer couplings J; nt = 30, 75, 160. The exchange interaction 
and anisotropy constants are Jfm = 100, Jafm = —50, D°^ FM = —1.0, D FM = 1.0. 
(b) Mean field theory (MFT) with the same parameters. 

ponents of each sublattice, thus allowing noncollinear magnetizations in both the 
FM and AFM layers. Our computer code, when specialized to a single magnetiza- 
tion direction, reproduces the results of reference [HH|. Without interlayer coupling, 
the code also reproduces the results for the monolayer ferromagnet and monolayer 
antiferromagnet simultaneously. The choice of anisotropies supports the orthogo- 
nal arrangement of the magnetizations of the FM and AFM layers favoured by the 
exchange interaction alone. The interlayer coupling destroys the perpendicular ori- 



82 



entation of the ferromagnet (in z-direction) with respect to the antiferromagnet (in 
x-direction), even at temperature T = 0, as can be seen from figure 1231 In this figure, 
we show the sublattice magnetizations calculated with GFT for three interlayer cou- 
pling strengths. With a positive interlayer coupling, all sublattice magnetizations 
develop a positive ^-component, whereas the x-components of the two sublattice 
magnetizations in each layer oppose each other. With increasing temperature, all 
x-components decrease until they vanish at a common temperature T^ 6el , slightly 
above the Neel temperature of the uncoupled AFM. For T > Tj^ el all sublattice 
magnetizations point in the positive z-direction. The AFM layer assumes a ferro- 
magnetic arrangement and remains so until a common critical temperature Tc is 
reached, at which the magnetic order vanishes altogether. 
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Multilayers 



5-layer (1 0-sublattices) antiferromagnet 



0.5 




temperature 

Figure 24: Green's function theory (GFT): Sublattice magnetizations of a 5-layer 
(10 sublattices) antiferromagnetic as function of the temperature. Parameters: 
J n = -100 Ji 2 = -30 J 22 = -86 , 66 J 23 = -20 J 33 = -73,33 J 34 = -10 J 44 = 
_60 J 45 = -8.66 J 55 = -46 , 66 D% = -1.0 (i = 1, 5). 

The model is easily extended to ferromagnetic, antiferromagnetic and coupled 
ferromagnetic-antiferromagnetic multilayers with individual parameters for each layer 
It is only a question of computer time. As an example we show in Fig. E3]the re- 
sults of a 5-layer (10 sublattices) antiferromagnet, where each layer has a different 
exchange interaction strength so as not to clutter the diagram. 

The theory could possibly serve as a basis for studying the exchange bias effect, 
where it seems, however, to be necessary to include interface disorder [T§1 |HUJ |H] 
in some way, for instance by introducing more sublattices per layer with different 
magnetic arrangements. 
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4.4. Working in the rotated frame 

In our exposition of GF-theory (e.g. see section 4.2.2), the higher-order Green's 
functions are all decoupled in a similar way independently of whether they are con- 
structed from S~ , S + , or S z operators or from mixed products of these. This might 
be a general weak point in the decoupling procedure — indeed, there is evidence that 
this democratic approach ignores essential differences in the roles of these operators. 
In particular, GF-estimates of the internal energy and specific heat are not as reli- 
able as those for the magnetization and it appears that this might be traceable to 
an inferior decoupling of the Green's functions transverse to the z-direction, which 
we normally choose to be in the direction of the anisotropy. 

Some recent publications [IHJ EH] suggest that working in a rotated coordinate 
system may provide a way to correct this deficiency, especially when considering 
the field-induced reorientation of the magnetization of a ferromagnetic Heisenberg 
film. The idea is that since the decoupling procedure for the single-ion anisotropy 
appears to function better in the direction of the magnetization than in the trans- 
verse direction, it ought to be better to change first to a rotated coordinate system 
where the decoupling can be carried out in the direction of the magnetization only. 
The angle of the rotation is determined from the condition that the commutator of 
S z ' with the Hamiltonian vanish in the rotated frame: [S z ', H '] = 0, where the 
prime refers to the rotated frame. This procedure is remarkably successful i n 
calculating the magnetic reorientation of a ferromagnetic film as a function of the 
external magnetic field in the presence of a single-ion anisotropy, as can be shown 
by comparing with the Quantum Monte Carlo calculations of Ref. jlH]. Not only 
that: the requirement [S z ', H '] — leads to an equation-of-motion matrix having 
no null-space — an enormous simplification of the entire calculation! 

Because of the apparent advantages of this new approach, we dedicate an entire 
subsection to it. First, we show how to implement the procedure, applying it to 
the exact treatment of the single-ion anisotropy; then, we present some of our own 
results and those of others |4TH I82j: finally, we discuss the method, examining the 
assumptions and pointing out some difficulties. 

4.4.1. The ferromagnetic film with an exact treatment of the single-ion 
anisotropy 

In this section, we show how to implement the GF-theory in the rotated frame for 
a typical case: the field-induced spin reorientation transition for spin S > 1. We go 
beyond the treatment in EH| in that we treat the single-ion anisotropy exactly 
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Consider the Hamiltonian (J2(Jlj) with a field B = (Bq,0, Bq) but without the 
dipole-dipole interaction and the K^-term. As the external -Bg-field is increased 
from zero, the magnetization vector initially in the z-direction rotates by an angle 
9 in the xz-plane, so that it points in the ^'-direction of a new frame (x',y',z r ). 
As in Ref. we shall do the calculations in the primed system, in which the 

magnetization vector has the components (0, 0, (S z ')). The transformation between 
the frames is 



/ (S x ) \ 




( cos e o 


sin 9 \ 




( (S x '} \ 








1 







(sy '} 


(252) 


{ (s z ) J 




y — sin 9 


cos# J 




{ (s z ') ) 





Because (S x ') = (S y ') = in the rotated frame, one need only calculate (S z ') in 
order to find the components of the magnetization in the original frame, once the 
angle 9 is known. 

To get the angle 9, an approximation is introduced: we demand that the com- 
mutator of S z ' with the Hamiltonian in the rotated system vanish. This implies 
that the following Green's function is zero: 

Gtf = (([Sr,H'];Sj')) = 0. (253) 

Evaluating the commutator yields a relation between Green's functions Gp~ = 

((St SJ '}} and G z +'~ = (((2S? - 1)#; Sy}}, 

(B$ cos 9 - B z sin 9)G ±'~ - K 2 sin 9 cos 6G'p~ = 0, (254) 

which, after applying the spectral theorem, produces the equation defining the re- 
orientation angle in terms of the corresponding diagonal correlations: 

{B% cos 9-Bl sin 0)C~ + - K 2 sin 9 cos = . (255) 

This is a generalization of the angle condition given in Refs 130] that can 
be used for the exact treatment of the single-ion anisotropy instead of applying 
the Anderson-Callen decoupling. Note that, as used here, the condition on the 
commutator must be considered an approximation. In Refs HO] the condition 
is fulfilled automatically because of the use of the Andersen- Callen decoupling. In 
general, the condition does not hold, as will be shown later. 

Following Ref. we introduce another approximation that in general also 

does not hold: we neglect all GF's not containing an equal number of S~ ' and S + ' 
operators. 

After transforming the Hamiltonian to the primed system and making the above 
approximations, the following Green's functions are needed: 

Gtj — (($t \Sj )) ) 
Gf" + '~ = (((Sr) n -\2Sr-l)St';S.')). (256) 
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The single-ion anisotropy requires that spin S > 1. Thus, in order to treat films with 
S = 1,3/2,2,..., one needs the first Green's function and those for n = 1,2,3,.... 
To get the equations of motion, the exchange interaction terms are treated by a 
generalized Tyablikov (RPA)-decoupling in which products of spin operators with 
equal indices are retained 



m rst ';sr '}} ~ r>«s fe + '-,sj '» + (st ')m r-,sj '» • (257) 



Now in the rotated system, (S£ ') = 0; i.e. the second term vanishes. After applying 
the decoupling procedure and performing a Fourier transform to momentum space, 
one obtains the following set of equations of motion. 



uG + >~ = 2(S Z ') + (S z ')J(g - 7k )G+>- 

+ (B* sin^ + B z cos 0)G + '~ + K 2 (l - § sin 2 6)G Z+ '- , 
uoG z+ >- = (6((S Z ') 2 > - 2S(S + 1)) - |J 7 k(6((^ ') 2 > - 2S(S + 1))G + '- 

+Jq(S z ')G Z+ '- + (B% sin^ + B Z cos 9)G z+ ~ 



uG {z)2+ >- = 8({S Z ') 3 > + 3((S Z 'f) - (AS(S + 1) - 1)(S Z '} - S(S + 1) 
+</7k(|S(S + 1) + (2S(S + 1) - 1)(S Z ') - l((s z ') 2 ) - M(s z f)G+>- 

+Jq(S z ')G iz)2+ '~ + (B*sin9 + B z cos6)G (z)2+ >- 



ujG {z)3+ '- = 10((S Z ') 4 > + 8((S Z ') 3 > - (6S(S + 1) - 5)((S Z 'f) 
-{AS{S + 1) - 1)(S Z ') - S{S + 1) 

+Jj*(lS(S + 1) + (25(5 + 1) - i)<S* ') + (3S(S + 1) - l)((S z ') 2 > 

-4((5 Z 3 ) -5((5 2 ') 4 ))G + >~ 

+Jq(S z ')G {z)3+ '- + (BjJ sin 9 + 5 2 cos 9)G (z)3+ '- 

+K 2 (1 - | sin 2 9) (2G (2)4+ '- - G (z)3+ '-) . (258) 



Here a = 1 is the lattice constant for a square lattice, q = 4 the number of nearest 
neighbours, and 7k = 2(cos/c x + cosk y ). 

As they stand, the equations (|258j) do not form a closed system. This, however, 
can be achieved by using formulas derived in Ref. [H] that reduce products of spin 
operators by one order (!), allowing the expression of some higher-order Green's 
functions in terms of lower order ones: 



+K 2 {1 - §sin 2 #)(2G (z)2+ '- - G z+ >-) , 



+K 2 (1 - | sin 2 9) (2G^ 3+ '~ - G^~) 



for S = 3/2: 



for S 



for S = 2: 



1 : 



G (,) 2 +,- + 3qz+,- ^ 



2 



2 



(259) 
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Inserting these relations into the system of equations (|258j) produces a closed system 
of equations. 

The equations of motion can be written in compact matrix notation 



(wl-T)G = A. (260) 



The quantities T, G, and A can be read off from equation (J258J), where the non- 
symmetric matrix T is a (2 x 2), (3 x 3), (4 x 4) -matrix for spins 5=1, 3/2, 2, 
respectively. The desired correlation vector corresponding to the Green's functions 

C = I ^ ^ I c?G1 5 

■ i {s- 'S^ n '\2S z ' -1)S + ') 1 ' 

is obtained via the spectral theorem. With the eigenvector method of Section 3.3, 
the components of the correlation vector C in configuration space are found to be 

d = f olkQik) = — dk x r alky Y, RifykLnAi (262) 
J n Jo Jo j,k,i=i 

(z=l,2,..,25), 

where the integration is over the first Brillouin zone and R(L) are matrices com- 
prising the columns (rows) of the right (left) eigenvectors of the matrix T and 
€jk = fijk/ie^ 3 — 1) is a diagonal matrix, in which uij are the eigenvalues (j = 1, .., 25) 
of the T-matrix. In sharp contrast to section 4.2.5, there are no zero eigenvalues of 

n 

Equation ()255|) and the set of integral equations ()262|) have to be iterated si- 
multaneously to self-consistency in order to obtain the magnetization (S z ') and 
its moments in the rotated system together with the reorientation angle 9. The 
curve-following method described in appendix B accomplishes this with alacrity as 
before. The components of the magnetizations in the coordinate system in which 
the magnetic reorientation is measured follow from the relations (|252J) . With the 
formulas from Ref. 07j it would be possible to treat the fourth-order anisotropy 
term — YU K^iiSfY exactly. A generalization to multilayers is also possible. 



4.4.2. Results of calculations in the rotated frame 

Here we describe results of calculations in the rotated frame, including results from 
the method described above. 

The paper deals with the Heisenberg ferromagnet with weak single-ion 
anisotropy in a varying transverse field. The Anderson-Callen decoupling is used 
in the rotated frame. The small anisotropies (e.g. for 5 = 2, K 2 = 0.01J) are 
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appropriate to 3d transition metals. For the reorientation as a function of the trans- 
verse field, there is excellent agreement with QMC calculations |48j. see Fig. EH In 
particular, the correlation (S x )/S is a linear function of B x , which is an improve- 
ment over calculations in the original coordinate system [2H], where the decoupling 
is performed for GF's corresponding to the components of the magnetization in the 
non-rotated frame. 

in rotated frame: A.C.-decoupling a ,QMC 9 
non rotated: (lines); S=2 J=1 00 K 2 =1 T=1 00 

1.0 
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Figure 25: Normalized magnetizations (S z )/S and {S x )/S and the reorientation 
angle for a spin S=2 Heisenberg monolayer for a weak anisotropy as a function 
of the external field: QMC (solid circles), Anderson-Callen decoupling in the 
rotated frame |49j (triangles) and in the non-rotated frame |28| (lines). 

If the anisotropy is treated exactly (see the previous subsection), the same thing 
is found, there being no difference from the results of Ref. for weak anisotropies 
within the line thickness. This astonishingly good result is perhaps the main point 
in favour of working in the rotated system. 

For the lanthanides, where values of the anisotropy can be of the order of the 
exchange interaction, the Andersen- Callen decoupling should break down and one 
would expect the exact treatment of the anisotropy to be superior. Surprisingly, the 
Anderson-Callen decoupling in the rotated frame still yields excellent results when 
compared with the exact treatment of K 2 and with QMC results [IB] for anisotropies 
up to K 2 < 0.2 J. This is seen in Fig. |2Hlfor the magnetic reorientation induced by 
the transverse B x -G.e\d for K 2 = 0.2 J and T = J = 100. The results of both Green's 
function theories (Anderson-Callen decoupled and exact treatment of the single-ion 
anisotropy) are nearly identical and deviate only slightly from the Quantum Monte 
Carlo results, which can be considered exact to within the statistical error. The 
reason for this is that at T = 100, the magnetizations from the two theories still lie 
very close to each other; at higher temperatures, this is no longer the case and the 
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results diverge beyond a certain value of B x . 



Compare: K 2 (exact) o , A.C.-decoupling A , QMC » 
S=2, J=1 00, K 2 =20 T=1 00 
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Figure 26: Normalized magnetizations (S z )/S and (S x )/S and the reorientation 
angle for a spin S=2 Heisenberg monolayer as a function of the external field: 
QMC (solid circles), Anderson- Callen decoupling 49J (triangles), present theory 
(open circles). 

Large differences must also appear as the anisotropy strength is increased, since 
it is known that the results from the Andersen- Callen decoupling do not approach 
the correct limit. This is evident from Fig. |23 where the field-induced reorien- 
tation for a Heisenberg monolayer with 5 = 2 from each GF theory is compared 
for a temperature somewhat below the reorientation temperature and for a large 
anisotropy K-2 = 0.5 J (T/J = 4.9). In this case, implementation of the Anderson- 
Callen decoupling along the lines of |4D] leads to a discontinuous transition from 
an angle (9 /(it/ 2) ~ 0.6) to full reorientation (9/(ir/2) = 1), whereas exact treat- 
ment of the anisotropy K2 produces a continuous reorientation transition. Such 
discontinuous transitions are also reported in Ref. |2U] fo r a treatment which is 
very similar to that of Ref. jlH]- The reason why discontinuities are not observed 
in Ref. j3H] is that only very small anisotropies are considered there. We at- 
tribute the discontinuous transition to the approximate Anderson-Callen decoupling, 
which is not justified for large anisotropy. The difference between the correspond- 
ing reorientation fields, Br, increases with anisotropy. For the present case it is: 

TV' exact , „ 

B^ 2 - Bft C - ~ 11 (for K 2 = 0.5 J). 

Unfortunately, we cannot say anything about the accuracy of the model treating 
the anisotropy exactly because there are no QMC calculations available for large 
anisotropies. The least understood approximation in this model is the generalised 
RPA decoupling of the exchange interaction terms of the higher-order GF, eqn. (|257|) . 
Previous calculations Jl] have shown (by comparing with QMC) that RPA is a good 
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Reorientation: K £ =50 T=490 J=1 00 S=2 




transverse field B* 



Figure 27: Normalized magnetizations (S z )/S and (S x )/S and the reorientation 
angle 0/(-7r/2) for a spin S=2 Heisenberg monolayer as function of a transverse field 
B x : Anderson-Callen decoupling 49 (dotted lines) and the present theory [oTi] (solid 
lines) for K 2 = 0.5 J (T/J = 4.9). 



approximation for a Heisenberg model (no anisotropy) with a field perpendicular 
to the film plane. To improve the present approach for a field in the transverse 
direction, one could resort to the procedure of which goes beyond the RPA with 
respect to the exchange interaction terms. 

We now consider a Heisenberg antiferromagnet monolayer with exchange anisotropy 
in a transverse field for S = 1/2. The Tyablikov decoupling is used. A recent paper 
[52"] reports results from an approximate GF treatment where the sublattices are 
rotated in such a way as to make the transverse component of the magnetization in 
each sublattice vanish. As in the ferromagnetic case, it is assumed that [5!f H'\ = 
at this angle, with the consequence that there are no zero eigenvalues of the result- 
ing equation-of-motion matrix. The authors describe their results as unexpected: 
the staggered magnetization of the easy axis shows a non-monotonic behaviour as a 
function of the transverse field and there is a nonvanishing easy-axis magnetization 
above the Neel temperature below a critical transverse field. 

To check the above results, we have computed the components of magnetization 
in the non-rotated frame directly from equations ()238|) of Section 4.3.2. Because we 
have developed [T3] a procedure to deal with zero eigenvalues of the equation-of- 
motion matrix, we do not need, contrary to Ref. |B2], any further approximation 
apart from the Tyablicov decoupling. In complete contrast to Ref. [82J, our results 
behave as one would expect: the easy axis magnetization decreases monotonically 
and vanishes as a function of the transverse field for temperatures above the Neel 
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temperature. Our results are shown in Fig. I2H We should welcome Quantum Monte 
Carlo calculations that could resolve the crass differences between these two sets of 
results. 

AFM monolayer in transversal field B z 

T= 0, 30, 45 T Neel (B z =0)= 47,2 




transverse field B z 

Figure 28: The magnetization components of a spin 5 = 1/2 antiferromagnetic 
monolayer (with the easy axis in x-direction) in a transverse field B z shown as a 
function of B z for different temperatures (T = 0,30,45). The Neel temperature 
TNeei(-B z = 0) = 47,2 for J=-100 and the exchange-anisotropy strength D x = —0.1. 



4.4.3. Discussion 

The most appealing aspect of decoupling in the rotated frame is not the excel- 
lent result for the field reorientation of the Heisenberg ferromagnet with single-ion 
anisotropy but rather the fact that the condition [Sf ', H '] — leads to an equation- 
of-motion matrix devoid of zero eigenvalues. Ref. may convey the impression 
that this condition is exact; if that were correct, decoupling in the rotated frame 
would undoubtedly be the method of choice because of the great simplification it 
offers. 

But "if" stands stiff. We offer a counter-example as a warning that the approx- 
imations used in Refs. [131301 |B2] should be taken with a grain of salt: an exactly 
solvable model demonstrates that [Sf ', H '] = is not in general valid! 

Consider a Hamiltonian having only an external field and a single-ion anisotropy: 

H=-Y, K 2 , k {Stf - Y,{B x Q St + B£S$. (263) 

k k 
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If it were true that a rotation angle could be found for which the commutator 
in the rotated system [Sf ', H '] vanishes, then the singular values of the matrix 
(SM'\[S Z ',H ']\SM) in the \SM) representation would also vanish. A numerical 
calculation for S = 1 shows that this is not the case for a finite K 2 . Furthermore, 
the numerical calculation also shows that the correlations (S~ 'S~ ') and (S~ ' S z ') 
do not vanish simultaneously with (S + ') = (S~ ') = 0. This shows that arguing 
with the Lehmann representation of the corresponding Green's functions as in Ref. 

is not correct because it is erroneously assumed that the intermediate energy 
states \m) are eigenstates of the z-component of the angular momentum. They are 
in fact, however, given by the superposition \m) = J2m c mM\SM), such that e.g. 
the relevant matrix element J2nm{ n \^~\ m )( m \^~\ n ) does not vanish in general. In 
Refs. |5(J| IH2] the Green's function G~-~ is taken into account correctly. 

Alternatively, consider finding a rotation that diagonalizes the model Hamilto- 
nian ([263)1 . If this were possible, then the commutator would be zero at the corre- 
sponding rotation angle, since two diagonal matrices commute. For this model, it 
is possible to show algebraically that no such angle can be found unless K 2 itself 
vanishes. 

At first sight it may seem strange that there is no angle at which the projection 
of the spin onto the z' axis is a good quantum number, for the non-commutativity 
of S z ' with H' implies that S z ' is not a constant of the motion but varies in time. 
But there is nothing wrong with this! One cannot argue that S z ' be time-invariant: 
the intrinsic anisotropy and the applied external field favour different directions and 
they do so according to completely different mechanisms. It would be wrong to think 
that there should be a "resultant" direction along which S z ' is quantized. Rather, 
the time-dependence of the opertor S z ' is simply a property of the Hamiltonian 
that must be respected. 

In conclusion, we regard the procedure of working in the rotated frame as not yet 
settled. It may in fact be advantageous if it succeeds in providing a more uniform way 
of treating the decoupling. The practice of employing [H, S z ] = is very likely much 
too severe in general. The method seems to work for the spin reorientation problem 
for the ferromagnet but is questionable for the antiferromagnet in a transverse field. 
The embedded null-space arising from a non-vanishing commutator is more likely 
an essential ingredient intimately bound up with the properties of spin. As such, it 
could be dangerous to ignore it. 
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5. Beyond RPA 

Up till now, with the exception of Section 4.2.5, we did not go beyond the Tyablikov 
(RPA) decoupling. In this section, we develop a formalism for treating the field 
induced reorientation of the magnetization for a spin 1/2 Heisenberg monolayer with 
an exchange anisotropy and, specializing to the magnetization in one direction, we 
show how higher-order GF theories discussed in the literature [SSJ E3 EHj follow 
quite naturally as limiting cases of our formalism. 

5.1. Field-induced reorientation of the magnetization of a Heisenberg 
monolayer 

We consider here a spin 5 = 1/2 Heisenberg monolayer with exchange anisotropy 
in an external field. We go beyond the Tyablikov (RPA) treatment by decoupling 
terms due to higher-order GF's. In the limit of the magnetization in one direction, 
we recover the results of Ref. for a vanishing anisotropy and of Ref. jHZj for the 
one-dimensional chain in the limit of a vanishing magnetic field. Without field and 
anisotropy one recovers the theory of Ref. [89|. The Hamiltonian under investigation 
is 



H — 2 51 ^ lm ^m^'t + Sm^l) 2 51 DlmS^Si ^Zi^ + $m + ^m) • 

Im Im m 

(264) 

The exchange interaction strength is Ji m , the strength of the exchange anisotropy 
is Di m and B ± = B x ± iB y with the external magnetic field B = (B x , B y , B z ). 

To get the equations of motion for the spin reorientation problem, the following 
first and second-order Green's functions are needed: 



^-ya-(l) _ I / QCi. Q—\\ 

G% (2) = (([£?, H];Sy)), 

GT = (([SfMSj)). (265) 



a 



The corresponding 8 equations of motion are 















{a = 4 


->->z 


^Gf x) 








uGf 2) 









(266) 
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The double-commutator Green's functions must be decoupled in order to obtain a 
closed system of equations. After Fourier transformation to momentum space these 
are 

(w-r)G q = I q , (267) 

a form that is amenable to the eigenvector method of Section 3.3. 

Generalizing the procedure of Ref. [HI] to the case where one has components 
of the magnetization in all directions of space, products of three spin operators are 
decoupled in the following way: 



SffiS? « a + ~cl*St + a z -cl+St + a z - C z +S z k , 

S?S+Si « afcjrs; + a + '<frS} + a-ttfSf , (268) 

where the correlation functions are defined as off = (SfSj). Here we have intro- 
duced the vertex parameters a + ~ , a , a z ~ and a zz , a +z , a~ z , where the indices 
refer to the indices of their associated Green's functions after the decoupling. In the 
limiting cases discussed later, we deal only with the magnetization in one direction, 
where only the vertex parameters a + ~ and a zz play a role. We show later how 
they may be determined by additional constraints. For the reorientation problem 
all 6 vertex parameters could play a role; however, for simplicity we assume that 
a + ~ ~ a ~ a z ~ and a zz ~ a~ z ~ a +z , in order not to have too many additional 
parameters. 

The inhomogeneities in eqn (J267|) are defined as the Fourier transformed ther- 
modynamic expectation values of the following commutators 



I(i = FT 







( \ 
q 


([Sf,Sf\) 




JZZ{1) 

q 


([Si,Sr]) 




/—(i) 
q 


([S z , Sj ]) 




p-(i) 
q 


([[St,H],Sr]) 




J+-(2) 

q 


([{S Z ,H},S Z }) 




jzz{2) 

q 


([[Sr,H],S.]) 




J (2) 


{ m,H],Sr}) ) 




{ k- {2) ) 
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2(S Z ) 



-(s-) 

zJ(l - 7q )(c+ - + 2c-) + zD{2c\l - 7q c+-) + 2^(,S 2 > + £+<S-) 
zJc+-(l - 7q ) + \B+{S-) + \B-{S+) 
-zJ(l - 7 q )cr ~ + z-DTqCio" - B-{S-} 
-zJ{l- lq )cl, -B-{S Z ) 



(269) 



Here, 



7q 



cos g for the linear chain with nearest neighbours z = 2 
| (cos + cos q y ) for the square lattice with z — 4. 

The T-matrix has the following form 



(270) 
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r 77 
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1^83 
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^87 
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(271) 



Without loss of generality, the external field may be chosen such that the reorien- 
tation of the magnetization takes place in the rcz-plane: B = (B x ,0,B z ). Then, 
(S y ) = and B + = B~ = B x , implying a number of symmetry relations for the 
correlation functions, such as (S + ) = (S~) = (S x ), cf+ = c[ m , cf z = c z ~ , etc. The 
non-zero matrix elements are then 



Tl5 

r 5 i 



53 



26 



37 



48 — 



■\B X B X - B Z B Z - \DB x z 1(l (S x ) 



+ zD\\ + a^{c z 2 l + {z -2)cg)) 

+zJD{\{\ - 7q ) + a + ~[(2 - 7q )(c» + (z- 2)cg) - \(d£ + (z - 2)c+f) 7q ] 

- a +-[(z-lh ?C --i(^-l)c+-]) 

+f(l- 7q )J 2 (l + a + -[2c- + c+- 

+(z - 2)(2c z 1 z 1 + c+f) - (1 + ^ 7q )(2c- + 4f)]) 



\B X B X + \DB x z lci {S x ) 



+ la + -Dj((c 20 - + (z- 2)c n -) 7q - { Z1 \ - l)c 
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10 



-f(l-7q)^ 2 « + (c 20 +(z-2)c u -c 10 (l + ^7q)) 
r 54 = 5*5* - HFz(l + 7q ) (S z ) + ^ 2 a + '2(^ - l)c^ 7q 

W" J/j([2 - z(l + 7q ) + 3 7q (^ - l)]cfo " c£ - (* - 2)c^) 

-z(l - 7q)</ 2 a + ~(c2 ~ + (* - 2)c^ - cfo (1 + *7q)) 
r 55 = 25 2 
r 58 = -25" 

r 62 = £W + zB x D(S x ) 1(l - zJDa zz 4 (l - 7q )(^7 q + 1) 

+§ J 2 (l - 7q ) (l + 2a" (c^ + (z- 2)4,- - (1 + *7q)cJr)) 
r 64 = -B^-B 2 + zDB x (S z ) + ^J J Da 22 (l - 7q )((z - l)c 2 ~ - (c 2 " + (z - 2)c 2 r)) 

-^J 2 (l - 7 q)« 2Z (^o + (* " 2)cfr - (1 + Z 7 qKo ) 
r 71 = \B X B X + ^DB x (S x ) lci 

+ i 2 DJa + -((c w - + (z- 2)crD7q - (^7 q - IKcT) 

-f (1 - 7q)^ 2 a +_ (c 2( T + (* - 2)c7f - (^ 7q + l)cr ") 
r 73 = -I^^ - B Z B Z - IDB x (s x ) 1(l + zD 2 [\ + a + -{c z 2 z + (z - 2)c 22 )) 

+zJD(±(l - 7q ) + a 4 " [(2 - 7q )(c 22 + (* - 2)c 22 ) - (z - l) 7g c 22 

+ |(^7 q - l)4f - |(4 - + (* - 2)c+f) 7q ]) 
+f(l- 7q )J 2 (l + a + -[2c 22 + 4 - 
+(* - 2)(2c 22 + 4D - (1 + *7 q )(2cg + 44)]) 
r 74 = 5*5* - zDB x (S z ){l + 7q ) + zD 2 a + ~2{z - l) 7q c 2 - 

zJDa + -([2 - z(l + 7 2 ) + 3 7q (z - l)]c 2 - - (c 2 " + (z - 2)c z u )) 
-z(l - 7q ) J 2 « + ^(c 2 - + (z- 2)c z n - (z 7q + l)cfo ) 

r 77 = -25 2 
r 78 = 25^ 

r 81 = ±B X B Z + IDB X {S Z ) 

+ IJDa + -(l - 7q )((z - l)c 2 - c 2 - (* - 2)c-) 
-f J 2 (l - 7 > + -(c 2 + (z - 2)c 2 r - (1 + Z 7 qKo ) 

r 8 3 = r 8 i 

r 84 = B X B X + zB x D(S x ) 1(l - zJDa+-4 -(l - 7q )(z 7q + 1) 

+§J 2 (1 - 7q )(l + 2a+-(c+- + (z - 2)c+f - (1 + «7q)4T)) 

r 85 = -b x 

r 87 = B x . (272) 

We have no explicit calculations with the eigenvector method for the spin reorienta- 
tion problem but we show now that, when specialized to one magnetization direction 
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only, limiting cases of the above equations lead to results found in the literature. 
5.2. Limiting cases 

For a magnetization in only one direction (the z-direction), the equations of 
motion reduce to a four- dimensional problem in energy-momentum space: 



10J 



r)G r 



Ir 



(273) 



with 



and 



q 

nzz(\) 

q 

(7+-(2) 

V cr (2) J 



(274) 



/ /q + ^ (1) \ 
JZZ(1) 

q 



j+-(2) 
q 

V if 2) ) 



2(S Z ) 


2B Z (S Z ) + zJ(ct - + 2cfg)(l - 7q ) + zD(2cH - 7q 4r) 



V zJcfo (! - 7q) / 

(275) 

5.2.1. Ferromagnet in a magnetic field, no anisotropy 

In this case, D = 0, B + = B~ = 0, B z ^ 0, leading to the theory of reference 
with the 4x4 T-matrix which is, in the notation corresponding to eqn. (J271j) , 



f 


r 5 i 



o 1 o \ 
o o 1 

o r 55 o 



(276) 



V o r 62 o o J 



where now 



r 5 i 

— ^2 



55 



«) 2 = -(1 - 7q ) J 2 (l + 2a"[{z - 2)4,- + c+T - (1 + z 7q )4T]) ; 
—B Z B Z + (cj+") 2 , with 

f(l " 7q)^ 2 (l + 2« + -[(z - 2)(ie+- + eg) 

+(^4r + <%) - (i + ^7q)(^4r + C fg)] 

2S 2 , (277) 



where z and 7 q are defined in equation ()270|) for the linear chain and the square 
lattice respectively For the linear chain (z=2), there are 7 unknowns ((S z ),cf^ , 

C 20 > L 10? ^20' 



cfn,Con,a + , oc zz ) and, for the square lattice (z=4), two additional unknowns 
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In both cases, the relations 

4r = 1/2 -(s*) 

eg = 1/4 (278) 

are valid. 

The eigenvector method of Section 3.3 yields 6 equations for the chain 

)i , 



c 



C j0 



zz 

So 



+~(zz) 

11 



P7T 

Cjo~ = \ J dqcos(jq)(B£Ll) : 
c% = \ jT dqcos(jq){K£IA) 2 (279) 

with j = 0, 1,2. 

For the square lattice, 8 equations are obtained (j = 0, 1, 2): 

Jo dqx Jo ^ C0S (^)( R ^ LI ) 1 ' 

P7V f"K 

± dq x dq y cos(q x j)(R£LI) 2 , 
Jo Jo 

/*7r rn 

±2 I dq x dq y cos(q x + q y )(~R£LY) m ■ (280) 

These equations do not yet suffice to determine the unknowns because the vertex 
parameters enter implicitly. In both cases, the missing condition is supplied by an 
expression for the intrinsic energy: 

Ei = ^ = -F(4T + eg) • (281) 
In order to evaluate (H), eqn ()115j) for (S~[Sf,H}) is compared with the explicit 
evaluation of (S^[Sf, H}) to yield 

zJ B z 1 1 _ i , . 

= "T-T + 2iv?S x < 282 > 

,~ M + B. + lJ(l- 7q ) . _ _ . 

e^-1 V q q ' 
The Green's function G^~ = G\ is the first component of the Greens function vector, 
given by 

Gi=£i2ii-2^. (283) 

Using 

Gi (u + z5) - G x (u) - id) = -2m ^ R u 5(u - ^) (LI), (284) 

i 

in eqn (J283|) . performing the ^-integration and comparing with eqn (J281j) yields the 
additional equation needed to determine all unknowns: 

-B z {S z )-%J(ct - + c? ) (285) 

= - Z -l-^ + U ^ \ u i + B* + l J(l-^) 

8 2 2 (i /;W Jt e^" 1 
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The equations (|278j) together with (J279l28Uj) and (J285j) determine the unknowns, 
from which one obtains the magnetization (S z ), the intrinsic energy Ei, the suscep- 
tibility x — d{S z )/dB z and the specific heat cy = dEi/dT. 

The numerical results in Ref. jBH| (obtained not with the present method but 
with the standard spectral theorem) demonstrate that RPA is a rather good ap- 
proximation for the magnetization and the susceptibility but that it is inadequate 
when transverse correlations play a role, as is the case for the intrinsic energy and 
the specific heat. In this case, it is very important to go beyond RPA. 

5.2.2. Ferromagnet with no magnetic field and no exchange anisotropy 
In this case, D=0 and B = and, for a linear chain [z = 2,7 q = cosg), one 
obtains the model discussed by Kondo and Yamaji 89 . Because of the Mermin- 
Wagner theorem, the magnetization is (S z ) = 0, and, for S = 1/2, isotropy demands 
that 

= |<&- ( 286 ) 

Therefore, one needs either the GF's ((S z ; S z }} and ({[Sf, H}; S z }} or ((£+; Sj)) and 
{{[Sf, H]; ST)). The first choice reduces the problem to two dimensions: 





TZz(l) 

q 

jzz{2) 

g 



(287) 



with 



Tf 



/-(D=0, Jf( 2 )=4c+-(l-cosg) , 
(1 - cosg)J 2 (l + 2a**(c&~ - (1 + 2cosg)cJ" ")). (288) 

These equations yield 

C zz ^ = ta. = h ( t ) (289) 

q u 2 - (u z q z{2) ) 2 2\uj zz \^uj - UJ ZZ U} + U ZZ J 

The standard spectral theorem produces 3 equations for determining the 3 unknowns 
aZZ i c io7> c 2o"- F° r spin 1/2, Cqq = (SqSq) = \. This, together with 

MT J 22 (2) Q 

c Z no = Ho = ~ I dqcos(nq) ^_ coth(||a;f|); n = 0,1,2 , (290) 

determines the unknowns. 

It is instructive to apply the eigenvector method of Section 3 and to obtain the 
same expression from 

Co = I ^ 7r dgcos(n g )(RfLl) i ; n = 0,1,2. (291) 

Here the matrix R consists of the right eigenvectors as columns of the nonsymmetric 
matrix in eqn ()287j) and L = R _1 consists of the left eigenvectors as rows and can be 
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calculated as the inverse of R. Note that R and L are not separately orthonormal. 
It is only necessary that RL = 1. One finds 




< 1 



2wf \ u * g -1 



/ 

8 = 









i 



(292) 



Evaluation of eqn (|291|) with these expressions produces eqn (|29U|) . 
Solution of these equations yields the intrinsic energy per particle 



^ ~ 2N ^ JnmjCnm + C nm) ~ 2^ C W ' (293) 



the specific heat per particle, 

dT 2 P d{3 



cy = — = \Jp—cfc, (294) 



and the susceptibility 



2 



X = PE C n=lE C nO- (295) 
n n=0 

The results of Kondo and Yamaji here reproduced numerically by the eigen- 
vector method, are largely in agreement with the exact calculations of Bonner and 
Fisher for a finite number of spins . 

The Kondo- Yamaji decoupling is generalized in Refs. [HHl ES] in order to treat 
the spin S=l antiferromagnetic Heisenberg chain. It is also used in Ref. j^Sl for 
the spin S=l low-dimensional quantum XY ferromagnet and in Refs. |94[ 19*5] for a 
kagome antiferromagnet. 

5.2.3. Ferromagnet with exchange anisotropy but no magnetic field 

This case (D 7^ 0, B = 0) is discussed in Ref. |HZj for the easy-plane XXZ 
chain, where 4 vertex parameters are used: af~ , «2~~, of, a^ z - These are fixed by 
the exact expression for the ground state energy and by assuming that the ratios of 
corresponding parameters do not vary with the temperature. 

A 2 x 2-problem results with the equations of motion 





r+-(2) 
q 

TZz(2) 

q 



(296) 



where 



and 



51 



J+-( 2 ) = 2J( C +- + 2 C -)(l- 7g )+2D(2c-- 7?C +-) , 
C (2) = 2Jc+-(l- Tg ), (297) 



2D\\ + at-c%) + 2JD(\{l- lq 
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+o+"[(2 - 7 „)cS - ic+- 7g ] + a+-[- 7g cg + K 2 7 g 2 - l)4f]) 
+ (1 - 7fl ) J 2 (l + a+-(2cg + c+T) - aj~(l + 27„)(2cS + 4")) > 
T 62 = -2J J Dc+ -af(l- T9 )(2 7? +l) 

+J 2 (1 - 7 ,)(l + 2of c+T - 2<(1 + 2 7g )c+-) . (298) 

The thermodynamics of the S > 1 ferromagnetic Heisenberg chain with uni- 
axial single-ion anisotropy using second-order GF's is treated in Ref. [HE]- The 
antiferromagnetic easy-plane XXZ-model for S = 1/2 is treated in Ref. [5%] . 
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5.3. The Tserkovnikov formulation of the GF theory 

Until now we have not considered the damping of spin waves. This is because 
we have neglected the influence of the self-energy, the imaginary part of which 
leads to damping effects. In this Section, we present a formalism which allows the 
treatment of damping. In the first subsection we develop the general formalism and, 
in the second subsection, we specialize it to a Heisenberg monolayer in an external 
field, evaluating the self-energy approximately. The formalism follows Tserkovnikov 
|96| . who derives a closed expression for the self-energy without making decoupling 
assumptions. For a review of the formalism, see e.g. Ref. [HZ]; a compact derivation 
can be found in the appendix of Ref. jHHI- The formal derivation of a Dyson equation 
for a Heisenberg ferromagnet is given in Ref. 

5.3.1. The general formalism 

The equation of motion for the single-particle Green's function G\ = ((A + ; A)) u is 

uG x = h + G 2 , with h = ([A + ,A]) and G 2 = (([A + , H]-, A)) . (299) 

Analogously, the equation for the two-particle Green's function G 2 may be written 

uoG 2 = h + G 3 , (300) 

with I 2 = ([[A + , H], A]). G3 is the three-particle GF, which Tserkovnikov ob- 
tains by a time derivation with respect to the second operator of G 2 as G3 = 
H];[A,H])). 

On the way to deriving an equation for the self-energy, Tserkovnikov introduces 
the ansatz 

(([A + , H\] A)) = C((A + ; A) + ((B; A)). (301) 
If one determines B such that ([B, A]) = the quantity C is determined by 

C = J a Jf \ (302) 

which can be proved by looking at ([[A + , H], A]) = C([A + ,A]) + ([B,A]). Intro- 
duction of the zero-order Green's function Go generates a generalized mean field 
expression: 

wG? = /i + / 2 /r 1 Go or G = -j—j—j ■ (303) 

A Dyson equation is now defined for G\. 

G x = G + G MG h (304) 
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where the mass operator M is defined by 

M = Ii 1 ((B;A))G^ 1 . (305) 



Gi = \ = —i , (306) 



The exact single-particle GF may then be written as 

h = h_ 

u - hli 1 - I\M io - I 2 I{ L - E 

where the self-energy is defined as E = I\M. 

The self-energy can now be expressed by 

E = (G 3 - G 2 G^G 2 )I{\ (307) 



The proof of this expression follows from eqns (jc 
E = uo - hG^ 1 - J 2 If 1 
= (G 2 GT 1 h-h)Ii 1 



G 2 Gl\G 1 u - G 2 ) + G 3 - ujG 2 )I^ 
= {G 3 -G 2 GfG 2 )I^. (308) 

5.3.2. The Heisenberg monolayer in an external field 

The Green's function G of eqn (|303|) leading to the generalized mean field expression 
for the Heisenberg monolayer in an external field (see the Hamiltonian of ([116)) ) is 
obtained from a Fourier transform to momentum space of 

h = ([S+,S;]) = 2(S z )5 ig 

h = ([[St,H],S;]) = 5 ig B2(S z )+ (309) 
+8 i9 Y: JiMsiSt) + (S+Sr)) -J2Ju5j2(S*S?) + (S+S g )] 



and reads 



with 



where 



2{S*)_ 
4 



G,, = ^ (310) 



K = B + i £( J q - Jk- q )(2^r + ), (311) 

c4^ (s ' s > ,q(R,_R,) ' (3i2) 



and 



TrY,((SD(S*} + ((s? - (SfiMSj - (s,))))^- 1 ^ 



2(^)(^)V + ^E^y eiq(Ri_Rj) - ( 313 ) 



N 
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With this expression, the dispersion relation is evaluated as 

E° = B + (S*)(J - J k ) + £(J q " Jk-q)^ + ) . (314) 

The first term corresponds to the Tyablikov (RPA) decoupling (see eqn (!126J) ). the 
term proportional to K* z corresponds to the fluctuations of the ^-component of the 
spin, and the term proportional to the transverse component is similar to the 
result of the Callen decoupling but with a different prefactor (see eqn (|141jl ). 

In order to describe the damping of magnons, one must go beyond this general- 
ized mean field approach, approximating the self-energy of eqn. (|3U8|) . The relevant 
term, which is the proper part in a diagram expansion j^Z] (leading to the name 
irreducible GF theory), is 

E(t) = Ga/r 1 = m,H]; [Sr,H)))^. (315) 

Evaluating the commutators yields 

MO = TTtL E JiiJ 9 A((StSl - SfSf); (S g S* - S-S*))). (316) 

A W ) lg 

A Fourier transformation to momentum space, together with the formulas of Section 
3.4 needed to derive the spectral theorem, allows one to express the self-energy in 
terms of the corresponding correlation function 



EH = — [°° J^—(e^ -1)1" dte iJt 

2% J-oo UJ — LO 



x E ■W„/ k:,r H ((-S .Sj - SjS* g ){StSt - St Sf)). (317) 



In order to proceed, the correlation function in this expression must be decoupled: 

/ / Q— QZ Q— QZ\f Q+ QZ Q+ QZ\\ 

* (s*s?)(s g -s+) - (s*s*)(s;s+) - (s z g sn(s;s+) + (s* g s?)(srs?) 

The longitudinal correlation function is approximated by its static value: 

(S-Sn(t)=^(t)=^(0). (319) 

The transverse correlation function is given via the spectral theorem by the single- 
particle Green's function G\ 

I roo ^—iu)t 

= 2^ J-oo dw e P u -\ \ i9 ^ U + ^ ~ Gl9 '^ U ~ i6 )) 

J_ 2(S Z )e~ lE ^ iq (R._R,) , v 

— N „ e^-1 • [ } 



q 
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Here, the full Green's functions G\ in the brackets have been approximated by the 
zero-order GF Go (a procedure which can be iterated to self-consistency) as 

2tlm-^- E G q , e- iq(Rl - R ^ = ^ E 2(S 2 >2tt<5(cu - ^g-^-R*) . (321) 

q q 

Now the t- and cj'-integrations in the expression for the self-energy (J317|) can be 
performed and, after a Fourier transform to momentum space, one obtains 

= 1 E - to ( J q - J^?K Z - (322) 

The single-particle GF is now specified and the magnetization can be calculated via 
the spectral theorem. The imaginary part of the self-energy describes the damping 
of magnons. This is the result obtained by Plakida in Ref. 

We are not aware of a numerical evaluation of the formulas above. The damping 
of magnons with the present formalism is, however, treated by analytical estimates 
for a two-dimensional S = 1/2 Heisenberg antiferromagnet in Ref. |lUlj and nu- 
merically in Ref. |lU(Jj and is also treated numerically in Ref. jHEl for a doped 
antiferromagnet within the t-J model. 
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6. Conclusions 

In this review we have given an overview of the formalism of many-body Green's 
function theory (GFT) and have applied it mainly to ferromagnetic, antiferromag- 
netic and coupled ferromagnetic-antiferromagnetic Heisenberg films. 

A prerequisite is that the systems to be examined have periodic structures in 
order to be amenable to the resulting two-dimensional Fourier transform from mo- 
mentum to configuration space. Any attempt to deal with local magnetic impurities 
would require calculations on a grid in real space, where one is limited technically 
by the number of lattice sites which can be taken into account. In this regard, 
the situation is the same as for Quantum Monte Carlo (QMC) calculations. A 
Green's function calculation including local magnetic impurities and a comparison 
with QMC results is reported in Ref. |lU5j . 

The crucial approximation in GFT is the decoupling of the higher-order GF's in 
the equation-of-motion hierarchy. The Tyablikov (RPA) decoupling yields reason- 
able results for the magnetization and susceptibility in one direction. This is seen 
by comparison of GFT with 'exact' QMC calculations for simple cases (see Sections 
4.1.5 and 4.2.5). If transverse correlations play a role, e.g. in calculations of the 
intrinsic energy or the specific heat, one has to go beyond RPA (see Section 5). In 
this case, third-order GF's have to be decoupled, requiring vertex parameters that 
have to be determined by additional constraints. For GF's of even higher order there 
is still no systematic procedure for the decoupling. Therefore, it is very difficult to 
make progress in this direction. In rare cases, e.g. for the single-ion anisotropy 
terms, it is possible to treat the corresponding terms exactly by using spin rela- 
tions that close the hierarchy of equations automatically with respect to these terms 
(Section 4.2.5). The exchange interaction and exchange anisotropy terms, however, 
have to be decoupled by generalised RPA procedures at the level of the higher-order 
GF's. 

A particular problem is the occurrence of exact zero eigenvalues of the equation- 
of-motion matrix. After application of the the spectral theorem, an adjunct term 
taking into account the corresponding null-space must be retained. If this term is 
momentum-independent, one can apply the standard spectral theorem in which the 
commutator and anticommutator GF's have to be used (Section 3.3). If, on the 
other hand, this term turns out to be momentum-dependent, the standard spectral 
theorem fails, and one must perform a singular value decomposition of the equation- 
of-motion matrix in order to eliminate the null-space from the matrix. This not only 
reduces the number of integral equations which have to be solved self-consistently 
but also makes the use of the anticommutator GF superfluous (Section 3.3). This 
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procedure is successful in a number of cases (see Sections 4.2.3 and 4.3.2). We were 
not, however, able to prove that this procedure works in general. For instance, we 
could not solve the spin reorientation problem of Section 4.2.5 in full because of 
numerical difficulties which we think are related to our inability to eliminate fully 
the momentum dependence of the terms connected with the null-space. We were 
able, however, to find an approximate solution for the spin reorientation problem 
with an exact treatment of the single-ion anisotropy by working in a rotated frame 
(see Section 4.4.2). The rotation angle is determined from the condition that S z ' 
commutes with the Hamiltonian in the rotated frame. This treatment simplifies the 
calculations because the null-space vanishes. However, we can show by a counter 
example that the condition above cannot be fulfilled in general. For the antiferro- 
magnet in a transverse field we can successfully deal with the null-space working in 
the non-rotated frame (see Section 4.4.2). 

The result of this full GF treatment deviates drastically from a GF calculation 
|H2~] in the rotated frame that employs the additional approximation that, as in the 
ferro magnet, the reorientation angle is determined by the condition [S z H '] ~ 
0. To resolve the discrepancy between these results, QMC calculations would be 
welcome. 

In most of the applications, we have considered only a simple square lattice 
because the double integrals in the Fourier transform from momentum space to 
real space can be transformed into a one-dimensional integral (Appendix C). This 
reduces the computer time considerably (by a factor of a few hundreds) because the 
Fourier transform must be calculated many times in the self-consistency procedure. 
There is, however, nothing preventing the use of double integrals directly for other 
lattice types if enough computer time is available. 

We are not aware of detailed numerical work which applies the Tserkovnikov for- 
mulation of GFT of Section 6.2 to Heisenberg films. This would allow a calculation 
of the damping of spin waves. 

We hope that we have succeeded in giving an overview of the present status of 
the application of many-body GFT to Heisenberg films that will stimulate the use 
of the reviewed techniques to related problems. 
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7.1. Appendix A: Calculating the intrinsic energy with GFT 

The following Heisenberg Hamiltonian for a monolayer is taken as an example: 

H = -B - 1 £ J«(Sr <?+ + S»5?) - 2 ^(^) 2 - (323) 
The intrinsic energy per lattice site is given by 

Ei = -B(Sf) -W M(S-S+) + (SfSi)) - K 2 ,((Srf). (324) 
z i 

In the following, we take S = 1/2 and S* = 1 as examples. 
1. S = 1/2 

In this case, the the single-ion anisotropy term is a constant because (S Z S Z ) = 1/4. 
In order to determine the intrinsic energy within GFT, one has to calculate the 
quantities entering eqn (|115|) of Section 3.7. Because = 1/4 and S^Sf = 

1/2 — Sf, the direct commutator yields 

Bp- = (sr[st,m~) = Be2-(s?))nAS!)-EM(s*s?)+(s-sf)+(s?srs+)). 

(325) 

A different expression for Bp~ can be obtained from 

St' = Jim £ s / ( G * + '~ (<* + *> " G -~ (<" " «)) ' < 326 > 

where G^' - is the Fourier transform of the Green's function 

Gtf = {{St;Sj)). (327) 

Equating the expressions (|325J) and (j326|) yields an expression for the intrinsic energy 
if relation (fSMj) . 

- £ J«<£J£?) = 2^ + 2S(5f> + £ 4<<^ + > + (328) 
is inserted into eqn (|325|) : 

Ei= ~ \B-{S!){\B + \zJ) + l -Y.MSlSiSt)-K,\ 

+ l\i m ly± f ^^(g+'-(lu + i5)-G+>-(lu-i5)). (329) 

k 

One now needs an approximation for calculating the Green's function and the ex- 
pectation values occuring in equation ()329|) . We have done this in Section 4.1.1 
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in the Tyablikov (RPA) approximation. The resulting Green's function is (see eqn 

r / \ (i S i ; Si ]) 2 ( S i) /oon\ 

G kl^J - -rpa ~ —rpa > K 6 ™) 

with the dispersion relation 

u>x pA = B + (S?)(J - J k ). (331) 

Now from eqn ()32fij) . 

+ ,_ _ 1 v f 2_{SQudw rpa\ _ 1 v 2(gf)^ 

The quantity J2i Ju(Sl $i St) lY1 equation ()329j) is obtained from the Green's function 

((S+;S*S-)), (333) 

which has same dispersion relation (|331j) but a different inhomogeneity ([St, S*Sf])\ 

1. e. 

((S + ; S z S-))^cu) = {{S ^ Sf ^} ) . (334) 

uj -uj k 

Applying the spectral theorem and a Fourier transform one obtains 

ET /qzq-QZ\ _ ^ 7 ( [jj ' Sj Sj ] ) _ 1 ^ T (jj ) fQQC^ 

JuVidi Si ) ~ e^ PA - 1 ~ iV k e^ PA - 1' 1 j 

We can now evaluate eqn ()329j) to obtain the following expression for the internal 
energy: 

Ei = ~\B - (St)(\B + \zJ) + i J>£™ - \J*) eP }5 A ) _ i - \K 2 , (336) 

which can be calculated after the magnetization has been determined self-consistently 
from equation ()130|) resulting from the spectral theorem in Section 4.1.1: 

(Sf) = ~ " (S^t) = ^-^E • (337) 

Knowledge of the intrinsic energy allows a determination of the specific heat and 
the free energy via eqns (jlllj) and ()112|) . 

2. S=l 

For 5=1 the single-ion anisotropy of eqn (j323J) is active. If the magnetization is in 
the ^-direction only, the exact treatment of the anisotropy of Section 4.2.5 requires 
the Green's functions: 

= ((St'i Sj )), 

Gtp- = «(2S? -l)St; S')) . (338) 
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with the exact equations of motion 

uGtf = 26 ij {St)-Y,M{{StSt-StSt)-,Sj)) + K 2 fitp-+BGtf, 

k 

ivGtp- = ^<65?5? - 4) + K 2ti Gtf -~£ ^*(«(65?5? - 4)5+; 5/}) 

+ 2((5 fc 5 + 5 + ;57))-2((^(^ + + 5 + ^;S7)))+5G^-. (339) 



We treat the single-ion anisotropy terms exactly, whereas we introduce RPA-like 
decouplings for the exchange interaction terms, taking care not to break terms with 
equal indices 

(((65f5f-4)5+;57>> ~ (65f 5f - 4>G+j- , 

((S^S^Sf;Sj~)) ~ (neglect of transverse correlations) . (340) 
A Fourier transform to momentum space yields 



uj-a - b \ GJ ~ \ _ / A+ - 
-c - d ) [ G z +>- ) ~ [ A z +>- 



(341) 



Here 



A+- = 2(S Z ) 

A z+ >~ = (6S Z S Z -A) 

a = B + (S Z )(J - J k ) 

b = K 2 

c = K 2 - l -{(SS z S z -A))A 

d = B + (S Z )J . (342) 

For a linear chain, J = 2 J, J k = 2cos/c; for a square lattice, Jo = 4 J, = 
2(cos k x + cos fcy). 

The eigenvalues of the matrix equations are 



1 

2" - ,v v 2 



w± = 5 + (5 2 >(J - |j k ) ± i/tff - -(6(5*5'*) - 4)if 2 J k + (-(S z ) J k ) 2 . (343) 



The Green's functions are then given by solving eqn (|341|) . 

+j _ _ A+<-(u-d) +bA z+ >- 
k ~ (u - lu+)(lu - lu-) ' 



Gk " (u-u+)(u-u-) ■ (344) 
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The spectral theorem then yields two equations determining (Sf) and (SfS?): 



(s-st) = 2 - (s?) - (s*s*) = 1 J2(ss + h = 4 E 1 



X 



A + '-(u + -d) + bA 



1 



; (A + >-(w- -d)+ bA z+ '-) — 



(S-{2St - l)St) = (St) - i(6<5?5?> - 4) = I £ ^ 1— ^ 



x 



{(^■>+ - a) + cA+n^J-j ~ (A^-(u- -a) + cA^)^^}- 

(345) 



Now, substitution of a) Ai = S { , Ci = Sf and b) Ai = S { , Ci = (2Sf — l)Sf into 
eqn f)l 15[) of Section 3.7 and insertion of the GF's (j344j) yields 



{(A+>-(^ -d) + kA+n-gj^ ~ (A + '-(-- ~d) + bA^n^hi) 



B 



z+, 



- a) + c^-)^^ - - a) + cA+O^i). 

(346) 

Calculating the commutators directly, inserting eqn (|324|) and eliminating J2k Jik({Si) 2 Sj:) 
by forming the difference 3B +, ~ — B z+ ~ and solving for Ei yields 



Et = |(35+ 
1 



B 



■--) - \{B + ff 2 ) - \{B + K 2 + zJ)(S?)) 



+ g E Jik \ ~ ( S k s f) - 2(S i Sf) + (S k (25? - 1)57 



(347) 



The first term comes from eqn ()346j) . Performing a Fourier transform on the last 
term gives 

H E M~3(S-S + ) k + (S~(2S Z - l)S+) k ). (348) 



This together with eqns ()345|) and ()342j) yields the final result for the intrinsic energy 
Ei = -\{B + K 2 ) - \(B + K 2 + zJ)(S z l ) 

+| E -T^{ [( 2 (^)(3(^ + -B- (St) J,) -K 2 + l(6S?S? - 4) J k 



+(6SfSf - A)(3K 2 - («+ - S - (5?>(J - Jk)] 



\ cu + - J k - 






'e^ + - l- 


+ 









}■ 

(349) 

For larger values of spin higher-order GF's are needed, but one can proceed 
analogously. The procedure applies of course to other Hamiltonians as well. 
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7.2. Appendix B: The curve-following procedure 

Consider a set of n coupled equations characterised by m parameters {Pf, i = 
1,2 ... , m} and n variables {Vf, i — 1,2, ... , n}: 

Si(P[m\;V[n\) = 0, for i = l,...,ra. (350) 

In our case, the parameters are the temperature, the magnetic field components, 
the dipole coupling strengths, the anisotropy strengths, etc; the variables are the 
spin-correlations. The coupled equations Si are obtained from the spectral theo- 
rem expressions for the correlations supplemented by the regularity conditions if 
necessary. 

For fixed parameters P, we look for solutions Si — at localised points, V[n], 
in the n-dimensional space. If now one of the parameters Pk is considered to be 
an additional variable V (e.g. the temperature), then the solutions to the coupled 
equations define curves in the (n + l)-dimensional space V[n + 1]. From here on, 

we denote the points in this space by {VJ; % — 0, 1, 2, , n}. The curve-following 

method is a procedure for generating these solution-curves point by point from a few 
closely-spaced points already on a curve; i.e. the method generates a new solution- 
point from the approximate direction of the curve in the vicinity of a new approximate 
point. This is done by an iterative procedure described below. If no points on the 
curve are known, then an approximate solution point and an approximate direction 
must be estimated before applying the iterative procedure to obtain the first point on 
the curve. A second point can then be obtained in the same fashion. If at least two 
solution-points are available, then the new approximate point can be extrapolated 
from them and the approximate direction can be taken as the tangent to the curve 
at the last point. 

The iterative procedure for finding a better point, V, from an approximate 
point, V°, is now described. One searches for the isolated solution-point in the 
n-dimensional subspace perpendicular to the approximate direction, which we char- 
acterise by a unit vector, u. The functions Si are expanded up to first order in the 
corrections about the approximate point, V°: 

Si(V) = Si(V°) + £ -j-AVj, (351) 

where AVj = Vj—V°. At the solution, the Si are all zero, whereas at the approximate 
point V° the functions have non-zero values, S°; hence, one must solve for the 
corrections AVj for which the left-hand side in the above equation is zero: 

£^A^ = -S°;{i = l,2,...,n}. (352) 
j=o uv i 
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These n equations are supplemented by the constraint requiring the correction to 
be perpendicular to the unit direction vector: 

n 

£«}AV} = 0. (353) 

3=0 

This improvement algorithm in the subspace is repeated until each of the S° is 
sufficiently small. In practice we required that J2i (S°) 2 < e, where we took e = 
1CT 16 . If there is no convergence, the extrapolation step-size used to obtain the 
original V° is halved, a new extrapolated point obtained and the improvement 
algorithm repeated. 

The curve-following method is quite general and can be applied to any coupled 
equations characterised by differentiable functions. By utilizing the information 
about the solution at neighbouring points, the method is able to find new solutions 
very efficiently, routinely converging after a few iterations once two starting points 
have been found. In addition, no single parameter or variable is singled out as "the" 
independent variable; instead, the (n+l)-dimensional curve can be viewed as being 
described parametrically in terms of the distance along the curve. This vantage 
point has great practical consequence: solutions in the neighbourhood of turning 
points (e.g. hysteresis for (S z ) as a function of field B) are just as easily determined 
as in any other region because the solution is always sought in a subspace nearly 
orthogonal to the solution curve. 



114 



7.3. Appendix C: Reducing a 2-dimensionsl to a 1-dimensional 
integral for a square lattice 

In the following we show how the double integral occuring from a two-dimensional 
Fourier transform when dealing with a square lattice (see e.g. eqn (|13Uj) ) can be 
transformed into a 1-dimensional integral. This transformation saves a lot of com- 
puter time in many of the applications discussed in the present review. 
Consider the evaluation of a double integral with the structure 



] /*7T /*7T 

I=— / /(cos k x + cos k v ) dk x dk v . 

TX l JO JO 

By substituting x = k x /n and y = k y /ir, this can be written as 

/ / /(cos nx + cos iry) dx dy. 
Jo Jo 



(354) 



(355) 



By making use of the fact that the integrand has the same value for all values of x 
and y satisfying the relation 



cos 7TX + cos iry = 27, 



(356) 



where 7 lies in the range (—1, 1), it is possible to reduce the double integral to a 
single integral over some suitable variable. The contours of constant 7 are shown in 
Fig. EH Each contour is given by an equation 



7iy(x) = arccos(27 — cos7rx). 



(357) 



Define now a function, ^(7), which is the area in the unit square in the xy— plane 



> 



"1 1 1 ( 
-0.2 -0.4 -0.6 -0.8 




Figure 29: Contours of constant 7 
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lying to the left of the line y(x) defined by Eq. 13571 for each value of 7. From the 
diagram, it is evident that the function ^(7) is given by 

1 f x o 

A{l) = — / arccos(27 — cos7ra;) dx (358) 

71 JO 

for 7 > and 

1 r 1 

A(j) = x\ H — / arccos(27 — cos7rx) dx (359) 

7T Jxi 

for 7 < 0, where x$ = - arccos(27 — 1) and x% = - arccos(27 + 1). These areas are 
shown in Fig. EH and Fig. |^ 

The double integral may now be written as a single integral over the variable A 
over the interval (0, 1): 

/= [ 1 f(2 1 (A))dA. (360) 
Jo 

In order to evaluate the integral numerically, it is only necessary to have an efficient 
representation of the function 7(A), so that a quadrature can be used to estimate 
the integral. A good strategy is to compute the function 7(A) at a sufficiently 
large number of points so that it can be accurately fitted to a cubic spline function. 
Thus, the labour involved in evaluating the integral I is enormously reduced, since 
the numerical representation of j(A) need only be computed once. 



> 




Figure 30: Area of the unit square to the left of a 7-contour for 7 > 



The numerical evaluation of 7(^4) is not without its problems, since the first 
derivative of the inverse function ^(7) has a singularity in its first derivative at 
7 = 0. Even though we need the inverse function, 7(A), whose derivatives go to 
zero at 7 = 0, there are still numerical difficulties in representing 7(A) by a spline 
function in the neighbourhood of A = 0.5; hence, it is better use a spline function 
to represent the function g(A), defined as (0.5 — A)j(A), and to get 7(A) from 
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> 




Figure 31: Area of the unit square to the left of a 7-contour for 7 < 

g(A)/(0.5 — A) in the neighbourhood of the singularity, using the value 7(0.5) = 
at the singularity itself. The function g(A), fitted to a cubic spline function, yields 
numerically stable values of the function and its first two derivatives. The spline fit 
to g(A) is obtained from values of the function tabulated at a set of knots equally 
spaced in the range (0, 0.5) plus values of the derivative of the function at A = 0, 
(derivative = — 1/(2-71")) and A = 0.5, (derivative=0). The second derivative of g{A) 
actually goes smoothly to zero at A = 0.5. Values of 7(A) in the range (0.5, 1) are 
obtained from the fitted values using the symmetry relation 7(0. 5+w) = — 7(0.5— u). 
The function 7(A) is shown in Fig. 1321 

While the above procedure allows one to obtain accurate values of 7(A) over the 
whole range of area, it does nothing to suppress the effects of the singularity in the 
first derivative of ^(7) at 7 = 0. These effects are not serious but they demand 
more effort from the integrator near A = 0.5. They can, however, be minimized by 
integrating out the constant part of the function f(2j), which is just /(0) itself: 

* = /(0)+ [\f(2 1 (A)-f(0))dA. (361) 
Jo 

We mention that in Ref. [HI] the double integral is transformed into an elliptic 
integral of the first kind with a transformation found in Ref. |102j . 
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Figure 32: 7 as a function of area A. 



7.4. Appendix D: Treatment of the magnetic dipole-dipole 
interaction 

In this appendix, we apply the generalized Tyablikov (|lfi5j) decoupling to the mag- 
netic dipole-dipole interaction. From this result the mean field approximation, 
as it is used e.g. in eqn (j!72|) . is obtained by neglecting the momentum de- 
pendence due to the lattice. After the decoupling procedure, the resulting term 
(([5?, H dipole }] (S*) m (£,~) n )) is added to the equations of motion (JUIJ), where H di P° le 
is the last term in eqn (jl60j) . After a Fourier transform to momentum space, one 
has the following additional terms in the equation of motion: 



' -n -ijt 

(T k y t+ 



—T£ \ ( G+' mn \ 
(T£)* G-' mn 
) \ G z r ; mn j 



(362) 



k 2± -ir^y 



where 




g(S z ) (T 20 + T ° 2 + |T 2 k + |T k 2 ) , 
fc<5*>(7*-7*+2i7*) , 




1Z± 

k 



^((^-)(r 2 k o + T k 2 -T 2 -T ° 2 ) 




(363) 



and 



T k = V 



I in 



X Z Urn 



exp(ik x xi) exp(ik y y m ) 



(364) 
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are oscillating lattice sums, which can be evaluated with Ewald summation tech- 
niques as outlined e.g. in Ref. 103 . 

This RPA treatment of the magnetic dipole coupling complicates the calculation 
of the magnetization considerably because of the presence of complex and dispersive 
(k-dependent) terms; therefore, we have neglected these terms in the applications 
and retained the non-dispersive terms only. This corresponds to a mean field treat- 
ment of the dipole coupling. In this approximation, eqns ()363|) reduce to 

T k = 0, 

n = f(5 + )(T 2 + T 2 ), 
T k 2± = -f(S->(T 2 ° + T ° 2 ). (365) 

This simplification takes the dipole coupling into account by an effective renormal- 
ization of the external magnetic field and leads to eqns ()172|) of Section 4.2.1. 

In order to justify this procedure we have done RPA calculations for the dipole 
interaction for two limiting perpendicular and an in-plane magnetization. 

In the appendix of Ref. [Hj, it is shown that, for these cases, a mean field calculation 
is a rather good approximation to the RPA result if the dipole coupling strength is 
much smaller than the strength of the exchange interaction, which is the case for 
many systems. We are not aware of a numerical treatment of the dipole coupling 
for the spin reorientation problem in GFT taking the dispersive and complex terms 
of eqn ()363|) into account. 

In the present review, we have applied the dipole-dipole interaction only in cases 
where the dipole coupling strength is small as compared to the strength of the 
exchange interaction, g/J << 1. Ref. jl 04) reviews dipolar effects in quasi-two- 
dimensional magnetic films, treating also cases g ~ J, g » J and J = with 
classical Monte Carlo simulations. 
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